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ABSTRACT 


This study investigates the optimal design of simple structures 
subjected to dynamic loads, with constraints on the structures' re- 
sponses. Previous studies have mainly dealt with static loads, and 
their methodology has been extended here to time dependent cases . The 
contributions of this work are in the formulation and satisfaction of 
the complicated dynamic constraints and in the insights gained into the 
nature of these problems. 

Three separate analyses search for the optimal design of ; (1) one- 

dimensional structures excited by harmonically oscillating loads, (2) 
similar structures excited by white noise and ( 3 ) a wing in the pre- 
sence of continuous atmospheric turbulence. The first problem has con- 
straints on the maximum allowable stress while the last two place bounds 
on the probability of failure of the structure. In all of these prob- 
lems, approximations are made in order to replace the time parameter 
with a frequency parameter. For the first, this involves the simple 
assumption that the steady state response is the area of interest. Xn 
the remaining cases, power spectral techniques are employed to find the 
root mean square values of the responses. The primary means of search 
for the optimal solutions is through the use of computer algorithms that 
combine finite element methods with optimization techniques based on 
mathematical programming. 
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A general conclusion is that the inertial loads for these dynamic 
problems can result in optimal structures t;.at are radically different 
from those obtained for structures loaded statically by forces of com- 
parable magnitude. In the case of the harmonically loaded structure, 
it is found that the design space can be disjoint. This makes the task 
of finding the global optimum difficult for even the simplest of prob- 
lems. 

An interesting feature of the optimal designs for cantilevered 
structures with a white noise excitation is that there is a tendency 
for some mass to be concentrated near the tip. The inertial forces 
from this mass tend to relieve the inboard stress. 

The wing in a turbulent gust environment demonstrates a possible 
practical application of the methods developed in the study. The model 
used contains a fuselage and nacelle and permits rigid body plunging as 
well as transverse bending. It is felt that the preliminary techniques 
developed are of practical value towards the design of aircraft that 
have fatigue life as an important design factor. 
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CHAPTER I 


INTRODUCTION 


A. PROBLEM MOTIVATION 

The goal of the field of structural optimization can be succinctly 
described as one of finding the structure of least weight that satis- 
fies certain specified constraints. The combination of more efficient 

algorithms with modern computers has expanded the capabilities of this 

$ 

field rapidly and to the extent that techniques have been developed 
that routinely optimize practical, statically loaded structures. Sim- 
ilar results for dynamically loaded structures have lagged behind due 
to the complications introduced by Inertial loads and the time param- 
eter. This thesis attempts a partial remedy of this situation by in- 
vestigating a series of dynamic response problems in order to find the 
least weight structure that can withstand the dynamic loads. 

The design of many engineering structures is influenced by the 
dynamic loads that act on the structure so that the search for opti- 
mal structures is a legitimate exercise. Landing impacts, gust ex- 
citation, rotating machinery and acoustic noise create loads on aero- 
nautical vehicles that are dynamic in nature and that are of primary 
importance in the design of aircraft substructures. Similarly, for 
as tronautical vehicles, rocket exhausts and atmospheric turbulence 
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are important design loads. These aerospace applications were the 
ones that were kept in mind while the methods of analysis used in this 
thesis were developed. However, other disciplines could benefit from 
the methods presented here. Specifically, for architectural struc- 
tures, earthquakes and winds create loads that are dynamic in nature, 
and these loads are playing an increasingly important role in building 
design. Further examples could be drawn from naval architecture and 
from mechanical design. 

Many of the examples mentioned above include loads that are sto- 
chastic, or random, in nature. Coupled with the fact that a large 
proportion of in-service failure of metal structures are due to fa- 
tigue, this provides a powerful motivation for studying the optimal 
design of structures under stochastic loading conditions. 

While no claim is made as to the direct applicability of the 
techniques developed in this work to practical problems in engineering 
design, techniques are developed and results achieved that could be a 
useful starting point for the more practical problems. The usefulness 
is enhanced by the use of constraints in the examples worked that are 
of practical interest in actual designs. For instance, constraints 
are placed on the maximum stress in the structure or on the fatigue 
life in the case of random loads. 

Due to the paucity of studies dealing with the optimisation of 
dynamically loaded structures, it is felt that this work makes sig- 
nificant contributions to the basic understanding of these types of 
problems. The inertial loads present in these problems can have an 


important effect on the loads a structure is required to withstand. 

The results obtained show that the optimal structure can be radically 
different from one obtained based on static strength considerations. 

B. RELATED WORK 

This section presents a survey of studies that have been done 
that relate to the present one,, pointing out their characteristics 
and how they compare with the present study. The thesis uses elements 
from a number of disciplines, but the unique portion of this work is 
the use of structural optimization and it is in this area that the 
survey will concern itself. Even in this specialized field, it would 
be impractical to give a comprehensive survey; instead, only the most 
relevant works are described. A more general view of the structural 
optimization field can be obtained from a survey article by Sheu and 
Frager (Ref. 1 ) while a text by Fox (Ref. 2 ) serves as an excellent 
introduction to the computational aspects of optimal design. Two re- 
cent conferences (Ref. 3 and Ref, provide "state of the art" de- 
scriptions of various portions of the field. 

Structural optimization with constraints on the dynamic behavior 
is a more specific area that includes the present study. A survey by 
Pierson (Ref. 3 ) on this subject divided it further into two subdivi- 
sions: (1) problems with eigenvalue constraints, and (2) problems 

with constraints on the response. The present investigation clearly 
falls into the latter category, but problems of the first kind played 
an important role in the development of the methodology used in this 
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report. In particular, references 6 through 10 are works that place 
constraints on the natural frequency or the flutter speed of the sys- 
tems to be optimized and theit provided a basis from which to attack 
the response problem. In fact, as Pierson pointed out, one of the 
dynamic response problems solved by Icerman (Ref, 11) has results 
identical to a problem with a natural frequency constraint that was 
first solved by Turner (Ref. 9)* 

It is to be understood that the five references cited above for 
the eigenvalue constraint studies are in no way inclusive of the con- 
tributions made to these problems. An attempt is made below to in- 
clude all the significant studies that have been conducted with con- 
straints on dynamic response quantities . These papers are a small, 
but rapidly increasing, part of the literature. 

The relative youth of the fieLd presents difficulties when one 
tries to classify the types of problems that have been studied. The 
ideal procedure would be to describe the problem that was studied, the 
method of solution and a discussion of results. Unfortunately, and 
unlike the more developed field of optimization with static loads, 
each paper treats a unique problem, usually in a unique way and, of 
course, obtains a unique result. Therefore, only the features of the 
studies that are relevant to the present work will be stressed in what 
follows . 

The youthfulness of the field is indicated by the fact that the 
earliest papers of this type known to the author were published in 
1968 . This work, published by Brach in two papers (Ref. 12 and Ref. 15 ), 
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found approximate optimal solutions for some one -dimensional struc- 
tures loaded by Impulsive or step forces. The problem formulation 
for these studies was in terms of finding the structure of fixed 
weight that minimized a specified deflection. This is a transforma- 
tion of the formulation used in this work: finding the structure of 

least weight with a constraint on the size of the maximum deflection. 

Fox and Kapoor (Ref, lU) published another "early" work that de- 
veloped a mathematical programming algorithm for finding the optimal 
design of truss-frame structures subjected to a half-cycle sine pulse. 
The response was estimated by using shock spectral techniques that 
gave conservative upper bound limits on the deflection and stress. A 
previous work by Fox and Kapoor (Ref. 15) made the important contri- 
bution of developing a simple technique for finding the derivatives of 
the eigenvalues and eigenvectors with respect to the system parameters. 

Levy and Wolf (Ref. 1 6 and RoL : . IT) provide a means of finding the 
fully stressed design for one-dimensional structures under dynamic 
loading. A fully stressed design is one where all structural elements 
exactly satisfy the stress constraints imposed on them. The motivation 
for their study comes from the fact that for determinate, statically 
loaded structures with constraints on the stress, the fully stressed 
design is optimal. For the impulsive loading conditions and the fi- 
nite element representations used, the solutions shown in these ref- 
erences are found to be optimal. However, fully stressed designs are 
usually not optimal in cases where more general dynamic loads are con- 
sidered. 
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A series of related papers by Venkayya, et al. (Ref. 18 and Ref, 
19) describes an optimality criterion that is used to find approximate 
optimal solutions for various types of dynamic loading. The criterion 
was developed specifically for problems with constraints on the natural 
frequency, so that it is exact for that case. When more general dy- 
namic conditions are considered, the results obtained have to be con- 
sidered as preliminary, qualitative designs. 

A specific area of practical Interest that can benefit from the 
methods of optimisation with dynamic constraints is that of the opti- 
mal design of structures to withstand earthquake loads. The 5th World 
Conference on Earthquake Engineering held in Rome in June, 1973 in- 
cluded four short papers on this topic. One of these, by Solnes and 
Holst (Ref. 20 ), replaced the dynamic load by an equivalent static 
load, so that it is not a dynamic response problem, strictly speaking. 
However, inertial effects are artificially included in the statically 
equivalent load. Another paper from the conference, by Nigam and 
Narayanan (Ref. 21 ), considered the excitation to be either a speci- 
fied acceleration or a probabilistic acceleration with a given power 
spectrum. The techniques employed in the paper and in another paper 
by Nigam (Ref, 22 ) to deal with the probabilistic nature of the ex- 
citation come closest to the techniques employed in Chapters IV and 
V of the present work to treat similar loadings. Another work of 
optimization for earthquake type of loads is given by Kato, et al, 

(Ref. 23). The loads in this case are approximated by shack spectra 
in a manner similar to that of Ref, lH. The diversity of models and 
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techniques used to study the optimization problem for civil engineer- 
ing structures indicates that it is a fertile ground for further re- 
search and systemlzation, 

A study that is more general in scope, but that has application 
to the earthquake problem, is contained in a recent report by Cassis 
(Ref. 24). In this case, the load is modelled as a half-cycle sine 
pulse and the response is obtained by performing a time integration 
of the equations of motion. The constraints considered include inte- 
grals of the time history of the response. This is one of the few 
papers dealing with dynamic response that retain the time parameter 
in an explicit form. It Is also the first report known to the author 
that includes mention of the fact that the feasible design space can 
be disjoint for certain types of dynamic excitations. This feature 
of such problems is one of the more exciting. The disjoint design 
space receives extensive treatment in Chapter III of this report. 

Chapter III deals with the optimization of structures excited by 
harmonically varying loads. In one sense, this is the simplest of 
the dynamic response problems since the time parameter can be removed 
by assuming that the steady state response is the only response of 
interest. By the use of energy methods, Icerman {Ref. 11) was able 
to develop an optimality criterion for one-dimensional structures ex- 
cited by a point load with an equality constraint on the displacement 
uxrectly under the load. In order to develop the optimality criterion, 
it was necessary to add the further constraint that the excitation fre- 
quency be less than the first natural frequency of the structure. Flaut 
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(Ref. 25) made a similar investigation but allowed the loading to be 
more general. While several examples were analyzed, and tlieir opti- 
mality criteria obtained, no explicit solutions were shown in this 
second study. Mros (Ref. 26) conducted a mathematically more rigor- 
ous study, which replaced the displacement constraints by one on the 
dynamic compliance of the structure. This is defined as the integral, 
over the entire structure, of the product of the magnitude of the load 
times the magnitude of the displacement under it. Despite the successes 
reported in these studies, the author knows of no effort that was made 
to expand on the results. An obvious, although difficult, extension 
would be to find an energy method that allowed the sinusoidal excita- 
tion to be applied at a frequency greater than the structures' first 
natural frequency . 

Finally, a series of papers that deal with static loads should be 
mentioned because of their relevance to the problem investigated in 
Chapters IV and V. They include some relatively early papers that 
sought optimal structures with constraints on their reliability (Ref. 

27 and Ref, 28). Moses and Kinser (Ref. 29) extended these results 
and used mathematical programming to find the optima. Araslanov (Ref. 
29) developed an optimality criterion that is applicable to simple 
beam structures loaded statically by forces whose properties are known 
only probabilistically. To do this, he defined the optimal structure 
to be the one where all cross sections have the same specified proba- 
bility of failure. These problems are the counterpart of the present 
study in that they assume that the structure and the load distributions 
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are described in some probabilistic manner but the loads are assumed 
independent of time. In the present work; the structural properties 
are assumed to be given and the loads ara constant in the space co- 
ordinate but vary in a probabilistic fashion with time, perhaps an 
enterprising investigator will integrate these two problems. 

C. SCOPE OF WORK 

The preceding literature survey omitted a few papers that were 
considered redundant or of little importance. It is quite likely that 
other papers were inadvertently overlooked. However; the survey at- 
tempted to demonstrate the full scope of the field of structural opti- 
mization with dynamic excitation and to indicate that this scope is 
still quite narrow. In addition; few of the papers cited were pub- 
lished; or, if published; were known to the author when this research 
began. For these reasons, the work reported on here does not build 
on the results of previous investigations to any major extent but 
rather attacks new problems. Of course; the tools needed for the 
analyses are gathered from existing disciplines, such as structural 
optimization, structural dynamics, aeroelasticity and probability. 

The core of the thesis is contained in three chapters that deal 
with three distinct optimization problems. In addition, a separate 
chapter describas the optimization algorithm used for the majority 
of the examples studied and an appendix details the finite element 
models that feed into each of the three problems. 
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The first problem is that of the structural optimization of one- 
dimensional structures excited by harmonically oscillating loads. 

This is similar to the cases dealt with in references 11, and 26, 

but a different approach is used that: provides added insight into the 
problem. In particular, the constraint that the first natural fre- 
quency of the optimal structure bo greater than the excitation fre- 
quency, which was an integral part of developing the optimality cri- 
terion of the previous studies, is removed. Another change is that 
the equality constraints on the displacements or the dynamic compli- 
ance are replaced by inequality constraints on the allowable stress 
within the structure. It is felt that these innovations provide for 
solutions of greater physical interest. Another facet of the present 
formulation is that the feasible region is disjoint. This provides an 
interesting theoretical result and one that may be of physical useful- 
ness as well. The major drawback c> t this formulation is that it is no 
longer possible to find an optimality criterion based on energy methods. 
This forces the investigator to deal with each problem on an ad_ hoc 
basis. One way of combatting this deficiency is the construction of 
analytical solutions to the optimisation problems by the use of con- 
cepts from optimal control. This rs a technique that met with some 
success when applied to problems with constraints only on the natural 
frequencies (Ref. 9 ). It is introduced, with limited success, in the 
present investigation because it holds the promise of providing solu- 
tions that analytically detail the effects of the various parameters. 
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Chapter IV deals with Jie second problem, which is the structural 
optimization of one-dlm£nsic*n.al structures excited by white noise uni- 
formly distributed along the span. A technical note by Nigam (Ref. Pi’ ) 
aided in developing the means for dealing with this type of probLetn, 
although the specific structures and constraints of Chapter IV differ 
substantially from those used by Nigam. Since the excitation is ex- 
pressed in probabilistic terms, the constraints also have to be evalu- 
ated using probability theory. Much of the chapter is therefore de- 
voted to defining the failure criteria used to evaluate a structure's 
lifetime. The methods ultimately used were obtained from Chapter 9 of 
a text by Lin (Ref, Jl) and include both fatigue failure and first ex- 
cursion failure. Further analysis in Chapter IV is devoted to out- 
lining how the response quantities and their derivatives, which are 
needed in the optimization procedure, are obtained through the use of 
superposition of natural modes. Finally, some numerical results are 
given and comments are made on points of interest. 

The methods of the earlier chapters are applied in Chapter V to 
a more practical problem, that of finding the optimal design for a 
wing excited by continuous atmospheric turbulence. The turbulence was 
represented by a power spectrum so that methods similar to those used 
in Chapter IV could be used to obtain the lifetime of the structure. 

A complicating factor is the translation of the atmospheric turbulence 
to the loads a wing experiences. A text by Bisplinghoff, et al . 

(Ref. 52, Chap. 5 ), provided the theory that permitted this. This 
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text also supplied the example (Ref. 52, Example 10.6) that was opti- 
mized, a tapered unswept wing that includes a nacelle and a fuselage 
and allows rigid body plunging in addition to wing bending. 

Finally, the last chapter summarizes the results obtained from 
the research and indicates areas that merit further study. 
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CHAPTER II 


OPTIMIZATION TECHNIQUES 

Tills chapter provides a description of the optimization methods 
that were used for the majority of the examples studied in this work. 
It does not attempt to describe alternative methods or to compare them 
with the methods used here. As mentioned in the introductory chapter, 
references 2 , 3 , and U collectively provide a good survey of the cur- 
rent state of various methods. 

Briefly stated, the methods used here involve coupling an inte- 
rior penalty function technique’ with a variable metric algorithm. 
These methods have been described elsewhere; in particular, Fox's text 
(Ref. 2) provides an excellent general presentation. Therefore, the 
present chapter gives only a brief outline of the method with emphasis 
on modifications developed during the use of the techniques. 

The first section defines terms that are common to optimization 
studies and are needed when the actual procedures are described in the 
following sections. A final section offers some observations on the 
algorithm based on experience gained from exercising it for the prob- 
lems of the thesis. 

A. CONCEPTS OF OPTIMIZATION 

The general process of optimization entails searching for the 
design that minimizes some specified function while satisfying all 
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the limitations applied to the design or its response. This section 
briefly outlines the concepts that put this general concept into quan- 
tifiable terms. Since finite elements are used for the majority of 
examples presented in this thesis, the development is presented in 
terms applicable to a finite element analysis. 

The first term to be defined is the objective or cost function. 
This is the function (or functional) to be minimized and is desig- 
nated by J . For the problems of this thesis, the cost function is 
always simply the sum of the design variables. 

The design variable is the second concept to be defined. This is 
an element of the system that may be changed in the process of seeking 
an optimum. The presenc study is concerned with one -dimensional thin 
walled structures whose design variables are the thicknesses of indi- 
vidual elements. The design variables are elements of a design vector 
that is notationally represented tv [t} .A related concept i.s that 
of the design space, which is simply the space of all physically pos- 
sible design variables. 

Limitations on the design arc termed constraints, and it is to 
the formulation and satisfaction of these constraints that the bulk 
of the effort of this work is directed. Tire constraints are desig- 
nated by the requirement that 

^ 0 , (i = 1 no. of constraints) } (2.1) 


- 14 - 


where is a function (exp licit or implicit) of the design vari- 

ables and the time and space coordinates. 

The simple two-dimensional example shown in Fig. 2.1 depicts 

these concepts plus some additional terms. This figure illustrates 

2 2 

the problem of minimizing the cost function J = t^ + tg subject to 
the constraints : 

8l « t x - 1 * 0 

S 2 = t 2 " °' 5 S 0 ’ 

g 3 = t l t a " 1 a 0 

The circular arcs are lines a Long which J is constant. The 
design vector is and the design space is given by all real 

values of and tg . This design space is divided into two re- 
gions by the constraint conditions: the "feasible" and the "infea- 

sible" region. The shaded, infeasible region is where the constraints 
are not satisfied, while the unshaded portion is the feasible region 
from which the optimal design must be found. While the optimal value 
of t^ = tg - 1.0 can almost be found by inspection (or by methods of 
ordinary calculus) in this case, it should be obvious that problems 
involving a large number of design variables and more complicated con- 
straints require considerable effort and ingenuity in the search for 
the optimum. 

k further concept that can be demonstrated with this two-dimen- 
sional example is that of the active and the inactive constraint. At 
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the optimum it is seen that constraints and g^ are satisfied 

as equality constraints (i.e., - g^ - 0.0 ). These are therefore 

designated active constraints. Constraint g^ is also satisfied, but 
the optimum does not lie on this constraint so that it is designated 
inactive . 

B. THE INTERIOR PENALTY FUNCTION 

When inequality constraints are imposed on the design, penalty 
function methods can be used to include the constraint in the objec- 
tive function. This strategem converts the problem to one that can 
utilize the powerful methods used to solve the unconstrained mini- 
mization problem. Reference £ contains a good description of these 
methods, and this presentation therefore focuses on the details of 
the particular penalty function used here. 

An interior penalty function is one that forces the trial design 
always to be in the feasible region. The specific function used in 
this work was of the form.: 

nc 

* = J - r 2 111 (S t ) • ( £ - 2 ) 

1=1 

The modified objective function, <& , is seen to become arbi- 

trarily large as the design vector approaches the constraint = 0 . 
As mentioned, this has the effect of forcing the trial design to be in 
the feasible region. Note that the form used here requires that the 
constraint be in the range 0 ^ g. £ 1 . This is accomplished by 
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redefining a given constraint so that it fits within these limits 
(e.g., the constraint t^ ~ 1.0 of Fig. 2.1 can be cast into 

the equivalent form g^ = 1.0 - . The r used in equation 

(2.2) is a specified scalar. The procedure followed is to minimize 
<J> for a chosen value of r and then to reduce r by some factor 
and repeat the optimization. In the limit as r 0 } the optimal 
result for the modified problem is seen to be arbitrarily close to 
the optimum of the unmodified problem. 

The extended interior penalty function method is a variation 
that was applied in reference 2h to a similar penalty function method. 
Figure 2.2 depicts an optimization problem that aids In explaining 
this refinement. In this diagram the function J - ax is being mini- 
mized subject to the constraint that x ^ b (or g - 1 - (b/x) ^ 0 ). 
While the modified cost function <I J - ax - r £n ( 1.0 - b/x) blows up 
as x approaches b , the extended penalty function remains finite. 
This is done by using the formulation: 

no 

® * J ' r £ G i( s i> > (2-3) 

i=l 

where : 


In (gjj 


g i * £ 


W - 


, S i ' £ 

In <: + g. < 6 
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The new expression is a TayLor series expansion of fn (g^) 
about the point in (g^) = e . 'll in reason for this esoteric con- 
struction is that the optimization process can now deal with designs 
in the infeasible region. Analyst:: of designs that are infeasible may 
sometimes be inadvertently performed either during the one-dimensional 
minimization described in the next section or by starting from an ini- 
tial design that is infeasible. 

The value used for e was selected by recourse to an argument 
similar to one used byCassis (Ref. £^) for a different penalty function. 
For the present penalty function, this gives e = exp {- r/$) . More 

comments on this choice for e are made at the end of this chapter. 

C. THE VARIABLE METRIC METHOD 

This section describes the particular mathematical programming 
algorithm used for the numerical optimizations of the thesis. When 
this study was in its early stages, a steepest descent algorithm was 
tried. The latter technique simply computes the gradient of the ob- 
jective function with respect to the design vector at the design point. 
A new design is then found by taking an improving step in the direc- 
tion of the gradient ("down the hill"). It is well known that the 
steepest descent method has very poor convergence properties but it 
was felt that it would suffice for the simple problems to be dealt 
with. For designs with more than three elements, this proved not to 
be the case. 


2C 



Following an investigation of various alternative algorithms, the 
variable metric method (also referred to as the Davidon-Fletcher -Powell 
method, after its developers) was settled upon. Reference 2 contains 
an excellent description of this method, and what follows is essen- 
tially a summary of that descr Lption . The variable metric method can 
be motivated by looking at a T'jylor series expansion, about design 
[Lq} , of an objective function: 

<£>(t) = $(t 0 ) + {>l(t 0 )} T {At} 

+ | {At] T pA(t Q )) {At} 

-f higher order terms , (2.4) 


where : 


n x l 

[At} = {t - t Q } 
n x l 

( a* ) 

f»(t 0 )] = — 

l3c i dt) - ft 0 } 

n x a 

[ a 2 > 

!^( 0 ) = 

l st i c ' e j Jft) = ft 0 ) 
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At the optimum [V<I>] a [oj , to that, to terms a C second order. 


near the optimum 


[v*(t )3 - {v*(t 0 )} + [v\(t c )] [At] - {0} . (P.5) 

Starting from { t^] , the indicated correction step is: 


[At] = - [^(t 0 )] _1 [VI>(t 0 )] 




If this were the actual procedure used to find the new design, it 

r J 

would be a second-order method. in practice, the [V 'i ] matrix is 
often difficult to obtain. Tills r.s particularly true lor problems 
dealt with here since the constraints used are very complex. 

The variable metric method was developed to circumvent this prob- 
lem by finding an approximation to the [V^O] ^ matrix. The method 
is outlined below and is followed by a brief justification of the al- 
gorithm. 

Directly from Ref. 2: 

(l) Start with some initial design vector { t] ^ and an initial 
positive definite matrix [HJ^ (typically the identity matrix) . Set 


(s } 0 = - [H] 0 (W(t 0 )} . 


(2) Find 

c (t , + a s ) 

q-l q q 


[At] = a^[s]^ , picking so as 

. The q subscript refers to the 


to minimize 


iteration number. 


(5) Compute: 


»Vi - w, + w, + - < s - 7 > 


where, designating 


[V} q = t^(t q+1 ) - W|, ' c q )5 

(A] q * r<»,(S},Cs]Jl/(l»]'{V} <1 ) 

{B} q - - lllH] q {v} q HlHl q {v) q } T ]/(Ev}^H] (i lv) q ) 


(U) Then set {s} . = ~ . r H‘! . ' V$(t .)} and return to (2) 

Cf + i. Q*ri CJ+i 


q+1 q 

until convergence is achieved. 


This rather complicated procedure can be heuristically justified 
by the fact that, for a quadratic, objective function with n design 
variables, the procedure yield:; [H]^ = [V^fr] ^ 

That is, if $ is of the form: 


<!> = {t} r tH3[l:} + [N ] [ t] + C 

Then: [HJ n == [M] 1 . (2.8) 


A proof of this statement can be found in Kef. 53* While the 
problems dealt with here are not quadratic in the design vector, the 
assumption is made that, close to tie optimum, the objective can be 
approximated by a quadratic, 
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1 . 


Interpolation. 


A remaining cask Is the eLajcnation of the rather innocent state- 
ment contained in step two of thi t lgorithm: "picking a so as to 

minimize 4(t + CL S ) 11 . This er tails performing a one-dimensional 

v q q q 

minimization at each iteration, am. it proved to be the most diffi- 
cult and time consuming aspect oE ihe optimization. The procedure 
finally settled on to perform this. 1-D search was a rather complex 

form of cubic interpolation that will be summarized here. 

* 


The goal is to find the value a that minimizes the scalar 
function ®(t + OS) . Assume that, the objective function can be 
approximated by a cubic equation in d 


<t» 


a h b3I t- (.Of 


^dO 5 


(2-9) 


If this approximation were axnct, the minimum could be readily 
found by setting the derivative of '1 with respect to OC equal to 
zero and then solving for d : 


O' o 

— = b + 2ca b 3dCT a 0 ; (2.10) 

cP 


* » 


* 


a 


c ± 



/3d 


. Ik - 


The choice or sign, is resolved by using the additional constraint 
that, at the minimum, the second derivative of the function is posi- 
tive : 


. 2 . 

b $ 

— - = 2c + 6d ot ^ o 

c V 


( 2 . 11 ) 


* 


Substituting the solution for a from equation (2.10) into 


(2.11) gives: 


6d 


2c + 


3d (- c± V' £ - 5db ') 


2 "1 
= ± 2 4fc - 5db > 0 


( 2 . 12 ) 


Clearly, the positive sign must be chosen. 

To complete the analysis, the values of the coefficients (a , 
b , c and d ) must be obtained. The original design is the value 
of the objective function at a ~ 0 . The slope of the objective in 

rn 

the a direction at a = 0 is given by {V' r, (t^)} [s] . Immediately 

then, a = ^(t Q ) and b = [ 1 C S } . The remaining coefficients 
are determined by evaluating 4 at two different values of Of . In 
order to assure convergence, these values were picked so that the min- 
imum was bracketed by the three function evaluations. 

Once Of is obtained using the above procedure, a test is made to 
see if it indeed is at the minimum value of S>(t+QS) . The test 
used was to compute j {V$ (t + Of S ) } T { s} j /( j {V$(t + a S))' 3 '[ | [s] J = Cl . 
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If a is exactly at the minimum. Cl is zero. The criterion used 
was that if Cl was less than some specified 5 then the optimiza- 
tion would proceed. If not, then an additional interpolation must be 
made utilizing the new values of <J>(t -fa's) and [V<l>(t +QI J 

until the criterion is satisfied. 

2, Minimum Thickness Constraints 

Under certain design conditions, it ia possible, in the absence 
of constraints on their size, that design variables may go to zero 
and even take on negative values. Since these design variables cor- 
respond to element thicknesses, it is physically and computationally 
undesirable for this to happen. Various methods have been constructed 
to deal with this problem, and this section describes a novel method 
used in Chapter III of this work. It is a method that worked quite 
well and is not well known in the structural optimization field. 

The technique used is a transformation employed by Pierson (Ref. 
10) for a continuous design variable. Modified to accept a discrete 
design vector, this transformation has the form: 

U3 - + \ {u S J . (2.12) 

The t , is a constant minimum thickness constraint while u 
mm 

is considered the new design variable. The beauty of this transfor- 
mation is the [t] remains positive even if {u} inadvertently has 
some negative components . 
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A minor difficulty arises when derivatives are needed with re- 
spect to the new design vector {u} . The indicated procedure is to 

use the chain rule by first taking the derivative with respect to [t] 
and then use 


&>] 

(a* 


( 

: - 


— ? 

N — 

Su ) 

Ut 

) 

( at 


( 2 . 13 ) 


D. COMMENTS 

Despite the analytical underpinnings described above, optimiza- 
tion techniques remain very much an art. It is felt that some per- 
sonal observation from one who began this work with a limited know- 
ledge of optimization techniques might prove of value to others who 
are in a similar situation. 

First, a disclaimer must be made to the effect that the use of 
the variable metric method coupled with an interior penalty function 
should not be considered a recommendation of either technique as the 
best method for solving a general problem. Each problem must be ap- 
proached on an individual basis, with a consideration of the require- 
ments and capabilities of each technique. The strong points of the 
method are that it is a sophisticated gradient method that proceeds 
to the optimum in a deliberate fashion. Other techniques, utilizing 
feasible directions (Ref, 3^0 or optimality criteria (Ref. 3)> are 
more efficient for certain applications and may even be better suited 
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for the problems worked here, A further general comment is that com- 
puter centers are now likely to contain optimization routines in their 
libraries. The first step for anyone embarking on an optimization 
problem should then be to determine if these readily available rou- 
tines are adequate or can be adapted for their needs. 

Given these general comments, specific perceptions gained while 
exercising the programs are offered below. 

The use of the - fn (g^ as a penalty function is an innovation 
with respect to structural optimization problems as far as the author 
knows. The more common interior penalty function is one of the form 
l/g , The log function seems to provide a smooth function with an 
easily calculated derivative. It would be interesting to hear of 
others' experience with different functions. 

The values chosen for the penalty parameter r of equation (2.2) 
have to be selected in an arbitrary manner. For this thesis, values 
of r ranging from ]0 to 10 ^ were used. The reduction - r^/lO 

was always used until the minimum r was reached. 

Texts on this method advocate iterating on each value of r until 
an optimum is reached before reducing it. This seems to be an unneces- 
sarily strict requirement and an alternative was used that reduced r 
after it appeared that little improvement would be made at the present 
value. This was done by specifying that if $ /$ - s l.o/(l.O + 10 r) , 

then r should be reduced by ten and the new optimization problem ini- 
tiated. If not, iteration continued until the criterion was met. This 
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approach has the added benefit that the criterion is satisfied quickly 
for large values of r and become}, increasingly more stringent as r 
is reduced. For the final value of r , a convergence criterion was 
employed. The criterion used was obtained from Ref, 2 and entailed 
checking if 


j{V$} T [H] {W}|/4 <0.02 

If the inequality xvas satis fied, the problem was considered 
solved; if not, the iteration continued. 

The use of the extended interior penalty function described in 
Section II. B proved to be of marginal value. The main reason for this 
is that the values of e were so small that the objective functions 
calculated using the extended penalty function were almost always too 
large to be of value in the interpolation procedure. This in turn was 
due to the way in which e is calculated. In order to assure that the 
transition point (i.e., e ) is between the minimum point and the in- 
feasible region it was found necessary to use & ~ exp (- r/$) 

Without going into detail, it is recommended that, if the extended 
penalty function is to be used, further efforts be made to obtain a 
better transition point when using the log penalty function, or that 
the l/g^ penalty function be used coupled with a transition point 
calculated by Cassis (Ref. 2h) : e = r/$ . 

A number of ether "tricks' 1 wore employed in the optimization and 
particularly in the one-dimensional search. However, it seems of 
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little value to detail them here. The main thing to be kept in mind 
is the nature of the optimization process and the mechanics involved. 
Some of the calculations of the next three chapters may appear exces- 
sive unless it is remembered why the optimization algorithm makes them 
necessary. 


- 50 - 



CHAPTER III 


HARMONIC EXCITATION 


A. INTRODUCTION 


Among the simplest dynamic response problems to formulate and 
solve are those of one -dimensional structures excited by harmonically 
oscillating loads. If only the steady state response is of interest, 
the time parameter can be removed from the equations of motion by as- 
suming that the structure responds at the frequency of excitation. It 
was supposed, therefore, that this type of problem would be a Logical 
beginning to a study of structural optimization in the presence of dy- 
namic loading. The results of this- chapter indicate that this suppos- 
ition is essentially correct but there there are unanticipated diffi- 
culties related to the fact that the feasible region is disjoint. In 
order to demonstrate this difficulty, some extremely simple examples 
are presented in the following paragraphs. 

Consider a uniform cantilevered rod excited by a uniformly dis- 
tributed sinusoidal torque. The differential equation and related 
boundary conditions for this system are (Ref. 32 , Chap. 3 ) 


S / se \ 3 2 8 

~ GJ 7 - ^ ~2 
a* \ ax / at 


- T e 


ift) t 
e 


( 3 . 1 ) 
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and 


£0 

9| X _Q S 0 atld GJ ““ 


= 0 


x=L 


Here, to is the excitation frequency. The amplitude of the 
steady state solution of Eq. (3*1) is 


0(x) = 

GJX 


cos yr x - 1 + tan L sin 


, (3.8) 


where 


X = 


2 

i or 

ce e 

GJ 


A graphical representation of 0(b) is presented in Figure $. 1 . 

The points to be made are that the magnitude of the deflection does not 
increase monotonically with the magnitude of the excitation frequency and 
that, given a specified deflection, there is not a unique value of the 
excitation frequency that results in that deflection. In fact, there are 
an infinite number of such excitation frequencies. This should provide 
an inkling of the problems to be encountered with a harmonic excitation. 
To make it more explicit, a further example Is presented below that in- 
volves a structural optimization problem with only two design variables. 
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B. 


TWO DESIGN VARIABLE EXAMPLE 


This section seeks the optimal design of a thin walled canti- 
levered rod excited uniformly in torsion by a harmonically varying 
load. If this system is modelled by two finite elements of equal 
length, equations from Section A.l can be specialized to the n = 2 
case to give the steady state equation of motion: 


2 2 
(CT L 
e OD 


24gj, 


0 


2(t 1 + t 2 ' 




(5.3) 


The constraints considered are that the magnitude of the stress 
be no greater than some specified value. From the Appendix, the stress 
can be expressed as: 



The motivation for representing the structure by two design vari- 
ables is that it is possible to depict the results graphically, thereby 
gaining a qualitative description of what would be encountered with a 
more realistic representation containing many elements. In this par- 
ticular case, an added benefit is that it is relatively easy to compute 


the stress amplitudes explicitly: 


2GR 


5T t„(l - X ) 
^ n 2 V e' 


max 


DET 


( 5 . 4 ) 


max 


2GR 

<®2 ‘ 9 1 ) “ 

L 


T (JX t* + t, (1 - 2X ) 
n w e 2 1 v e / 

DET 


vmere: 


2 2 
CD I JL 
e 00 

24gj„ 


t = tlr/(4j. s ) 

n ' ' 0 max' 


DET = - Sk e f - t a (6 > e - 3^)} 


The constraints for this problem are that the absolute values of 

S./S must be less than unity. In the notation of Chapter II, these 
1 U12.X 

are written as g. = 1 - (S./s s 0 

i v x max 

Figures .2 and 3*3 show the feasible and infeasible regions for 

values of X^ ranging from zero to three and for equal to 0 . 085 . 

For the X^ = 0 case (static loading), presented in Fig. 5*2(a), 

the constraints are seen to be two straight lines. The cost function 

is simply J = t^ + tg so that the optimum is at the intersection of 

these two lines. Figure 3*2(b) shows the design space for X = l/2h 

6 
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X 
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FIG. 3- 2- -Design Space for a Cantilevered Rod Excited 
in Torsion at Nondimens ional Frequency X . 






Til is is 


and it is seen that there are two separate feasible regions, 
the difficulty that is illustrated by this example and is discussed 
further below. 

Figure 3*3( a ) shows the results for X^ = l/6 , The constraint 

at t^ = 0.0^ is a minimum thickness constraint that is included to 
eliminate the t^ =0.0 solutions that satisfy the stress constraints 
but are physically unrealistic. It is seen that the least weight solu- 
tion is in the upper region at t^ = 0 ,0^ and tg - 0.23 . Finally, 

Fig. 3-3(b) shows that for the = 3 case there is again only one 
feasible region. 

The explanation for this curious behavior is to be found by study- 
ing the eigenvalues of the system. Let X^ and X^ denote the non- 
dimensional values of the first and second eigenvalues. In Fig. 3*3 ( a )j 
the designs with X^ equal to the excitation frequency are all on a 
straight line emanating from the origin with an equation given by 
t^ = 2.06 tg . This line proceeds directly through the middle of the 
infeasible region, dividing the design space into two distinct regions. 
Clearly designs that have X^ = X^ are infeasible because this repre- 
sents a resonance condition with an unbounded response. The region 
> 2 .06 tg contains designs where X^ > X^ while the region 

t., < 2.06 t_ contains designs where X < X . Each of these 

12 1 e 

regions has its own optimum, as is demonstrated by the figure. 

Segenreich and Rizzi (Ref. 35) have shown that the eigenvalues of 
cantilevered rods modelled in the fashion described in the Appendix 
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have prescribed limits. For the specific case of two design variables 
of interest here, these limits are given by: 


0 £ x x <; 0.5 , 

0.5 £ V, * 2.0 . ( 3 . 5 ) 

As a function of the excitation frequency, there are, therefore, 
either one or two distinct feasible regions. These regions are given 
by; 


Excitation 

No . of 


Eigenvalue 

Frequency 

Regions 


Relationships 

o 

ii 

0 ) 

r< 

1 


X- > X > X 
2 1 e 

o < x £ 0.5 

2 j 

i 1 *) 

0< \<\<\ 

e 


la.) 

0 £ \ < X e < X 2 

0.5 < X <2.0 


(i.) 

X 1 < X e < X 2 

2 



e 


(2.) 

X 1 < X 2 < X e 

X > 2.0 

D 

1 


X e > X 2 > X 1 

The X >2*0 case 

explains why 

Fig. 

3.3(b) contains only one 


feasible region; the excitation frequency is greater than any possible 
eigenvalue of the system. 
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It should be clear why this disjoint property of the feasible re- 
gion presents a difficult obstacle in the search for a global optimum* 
While it is possible to analyze the two design variable case graphi- 
cally in a thorough fashion, this is not practical for designs with a 
greater number of elements. Figure 3*1 was presented to motivate the 
hypothesis that for the continuous case there are an infinite number 
of local optima corresponding to the infinite number of distinct re- 
gions where ^ , i = 1,2, . . 

For problems with an arbitrary number of elements, some method 
such as that described in Chapter II has to be utilized to search for 
an optimum. But such methods have the drawback that the search takes 
place inside one feasible region. Therefore, for a given problem, the 
global optimum is found by selecting the minimum of all the distinct 
local minima. Cassis (Ref. 2U) encountered disjoint feasible design 
spaces while studying a different dynamic response problem and found 
it preferable to search for the optimum in the infeasible region by 
using an exterior penalty function method. His thought was that the 
solution would be more likely to proceed to the global optimum. But 
this technique provides no advantage here since an exterior penalty 
function technique still proceeds "downhill" and would not, therefore, 
cross over the infinitely high "ridge" where the excitation frequency 
equals an eigenf requency in order to descend into the "valley" of the 
global optimum. More comments are offered on this problem in Section D. 
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It might be supposed that the disjoint nature of the feasible re- 
gion is due to the omission of structural damping; in a sense,, this is 
true. The addition of damping gets rid of the infinitely high ridges, 
since a damped structure excited at its resonant frequency has a fi- 
nite response. A brief study that included damping vias made, and a 
result from the study is presented in Fig. 3 *^- The figure superim- 
poses the k “ l/6 case of Fig. 3 .3 ( a ) and the results from an iden- 
tical problem except that the shear modulus was multiplied by (1-fia) 
where a is a small structural damping factor. This is a technique 
frequently used to take account of the fact that structures have damp- 
ing present in them (Ref. 36, Chap. 12 ). The value of a used to ob- 
tain the results shown was 0.1 — an unrealistically high value, but 
one that depicts the damping effect clearly. It is seen that the 
damping reduces the infeasible region and prevents it from extending 
to infinity. The disjoint character of the design space has been 
eliminated, but two minima are still retained as pockets of the uni- 
fied design space. The basic problem of finding the global optimum 
still remains. Note that the optimal solutions for the damped case 
do not differ greatly from the undamped case. Damping was not in- 
cluded in the analyses presented in the remainder of this chapter 
since it was felt that the benefits gained from added practicality 
or realism do not offset the complications introduced by complex 


variables 




c. 


FUNCTION SPACE SOLUTIONS 


Before proceeding to the finite element solutions, another pro- 
cedure that is applicable to these sorts of problems is presented: 
that of solving optimization problems by dealing with the differen- 
tial equations directly. The motivation for this section comes from 
the success others achieved while applying optimal control techniques 
to structural optimization problems. In particular, Weisshaar (Ref. 7) 
and Armand and Vitte (Ref. 8) were able to find optimal thickness dis- 
tributions for a number of problems that had constraints on the system 
eigenvalues . 

This section develops the criteria for an optimal solution for a 
harmonically loaded structure and solves some special cases. 

Only one -dimensional structures are used in this study; therefore, 
the equations can be put into the first order form generally used in 
control theory: 


[*}' « U(t,s)]{x} + £ P} . (3.6) 


With boundary conditions at s = 0 and s = 1 . The terms used 

are defined as: 


[x] 


(x(s)} - n x 1 vector of state variables 


t = t(s) = thickness distribution, the control 
variable of the problem 
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and 


[P] = n X 1 vector of the load amplitude 

( )' = denotes a derivative with respect to s 

s = the nondimensional coordinate and inde- 
pendent variable. 

The analysis given below is an application of the methods de- 
scribed by Bryson and Ho (Ref. 37). Only the barest outline of the 
procedure is presented here. 

The problem statement used in this thesis is that of minimizing 
the weight subject to constraints on the response. Mathematically, 
minimize 

1 

j = J tds . (3.7) 

0 

Subject to 


{g(x,s,t)J £ o , 

q X 1 vector . (3.8) 


The Hamiltonian is constructed by using standard procedures of 
Ref. 37 


H = t + (\} T ([FlU] -i- {pj )+ {(Jl] T (g} , ( 3 . 9 ) 
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where 


{X} = n X 1 vector of adjoint states , 

[p] = q x 1 adjoint vector for the constraints 

The value of is zero when j 0 and is so when = 0 . 

The Euler-Lagrange equations are: 


( ] 

m t ( as 1 


M' - - - 

- - [x} T [f] - {p} T I — 

. (3.10) 

( a^ ) 

(ax J 



And the "control equation" is; 


The transversality condition provides the required boundary con" 
ditions : 


+ [p] 


SB 

at 


(3-11) 


an 

at 


= 0 = 1 + i>} 


a? 

at 


£*} 


M 


T 



(3.12) 


It is felt that the method is best dealt x-fith here by example. 
Hopefully, these examples also clarify the technique. 



1 . Example: Cantilevered Beam With a Static Load 


A cantilevered beam acted upon by a uniform static load has a dif 
ferential equation and associated boundary conditions given by 


d 


2 




= P 


(3.15) 


dw 

1 .2 
d w 


d 

/ d w \ 

^ — 

= El — 


— 

EI — 

dx 

n dx 

'x=0 

x=L 

dx 

\ d* 2 / 


x~L 


0 


Using notation and assumptions given in Ref. 7, the first order 
form of this system is 





where 


x 


1 


= w/l 


> 
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and 



1 

dw 

X 2 * 

L 

dx 


t 

,2 
d w 


1? 

7~2 

dx 


1 

d 

= 

? 

dx 

p - 

PL 3 

1 


,2 
d w 


dx 


El, 


The optimization problem is specified as that of finding the 
thickness distribution that minimizes the total weight while satis- 
fying the constraint that the magnitude of the stress along the span 

of the beam is less than some specified S . By use of the fa- 

2 2 

miliar formula S = (Ed/2)(d w /dx ) , this constraint can be put 

into the form: 


S 1 



( 5 . 1 ;) 


O 

where a = EdL / 2S 

' max 

Equation (3*9) is then: 


X 2 X 3 

H = t + 




!"2 ’ + V 1 * * V + 


(3.16) 
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Equations (3.10} and (5*11) are evaluated to give: 



cJH 

— = 0 
c>t 


X v, H x , 


_ 23 


i ^i + - 

t t 


(3.18) 


The last substitution is made because if |u ^ 0 , [a(x^)/t] = 1 

The notation sgn ( ) designates that only the algebraic sign 
of the quantity is used. 

The boundary conditions on the adjoint variables are 


\(i) = V>u> 5=1 V°> = V°> = 0 

The first order equations and the boundary conditions give immediately 

= ?(s - 1) , 

Xj = P(s 2 - 2s 4- l)/2 } 

\ = 0 ' 

= o . (3.19) 
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These results can be placed in Eq. (3*15) to give: 


1 + £ = 0 t - - p 


0 . 20 ) 


Since p is zero only when a|x^j/t is less than unity, it is 
clear that ^ (1, - 2s 4 s 2 ) . Equation 0-20) states that if 

p is zero, then t = 0 and the inequality on t is violated except 
at s = 1 . Therefore, p cannot equal zero, requiring ajx^j/t = 
1.0 across the span. Stated another way, this says that the optimal 
solution is the one that creates a fully stressed structure. This is 
a well known result for problems of this type with a static loading. 

The entire solution can now be written down as: 


t = - p = aP(l - 2s -t s 2 )/2 

x 2 = s/a ^3 = aS 

x x = s 2 /2a \ = as 2 /2 . (3.21) 

The ease with which this analytical solution was obtained makes 
it appear that solutions with a harmonically oscillating load might 
also be tractable. The formulation for the same problem as above 
exceDt that the excitation is harmonic with frequency cn can be 
written in terms of the static problem by adding several terms. The 
subscript ( ) in the following equations refers to the static case 

S L 
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\ 


and r is a nondimens ional frequency equal to ,n c W /l % * 

With this notation* the changes in Eqs . (3 .1^ )- (3 . IB) fur the harmon- 
ically excited structure are: 


H 


at 


r^tx. + p 

e 1 


h + Xi r tx. 

st 4 e L 


- r e tX l, 


au 




st , ,.2 

- + V e x l 


a 0 


0 


) 


These additions prevent determining that the optimal structure 
is fully stressed. Without this* it has been found impossible to 
treat these equations analytically. Numerical techniques that solve 
the two point boundary value problem and the associated control equa- 
tion have been applied with little success. The main difficulty is 
in dealing with the stress constraint. The character of the solution 
changes at the value of s where the constraint changes from being 
inactive to being active. This requires patching together arcs as 
explained in Ref. 37* Chap. 3* If numerical techniques are to be 
used, it seems preferable to convert the problem to an unconstrained 
one by using the penalty function method as described in Section II. B. 
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If this is done, Eq. (3.7) becomes 

1 
J 

0 

and the difficulty in patching arcs is avoided, albeit the formula- 
tion becomes slightly more complex. 

2. Example: Torsional Rod Excited by a Harmonically Varying End Load 

Consider a torsional rod that is being excited at its tip by a 
harmonically oscillating load with frequency CD and constant ampli- 
tude T . Pose the problem of finding the thickness distribution 
that minimizes the weight of the structure subject to the constraint 
that the tip rotational amplitude is equal to a specified value D 
This problem was first solved by Iceman using energy considerations 
and with the additional constraint that the first natural frequency 
of the structure be greater than the excitation frequency. It is 
similar to a problem studied by Ashley and McIntosh (Ref. 6) and by 
Turner (Ref. p) who found the minimum weight structure for a canti- 
levered torsional rod with a fixed tip mass and an equality constraint 
on the first natural frequency. 

With the familiar assumptions that GJ = GJ Q t and I Q = , 

the differential equations can be put into the form (Ref. 8): 


/ 


[t - r in (g)] ds 



(3.23) 
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And the associated boundary conditions are 


where 


i x 1 (0) - 0 

x^(l) - D 

x 2 (l) = T 


r 2 a 


P 

a> f T L 
e QO 

GJ„ 


2 


Note that the equality constraint and the excitation are con- 
tained in the boundary conditions. The inertial loads appear in the 
_2 

- r t term. 

The Hamiltonian of Eq. (3.9) has the form: 


^1 X 2 2 

H = t •!• - XgP tx x 


(3*2 ! 0 


Equations (3.10) and (3. 11) give the relations 


0 r c t' 


- i/t 0 


(5.25) 




— « 0 * 1 - 

at t 


X 1 X 2 , „2 


2 - X 2 Fx 1 


(3.26) 


and the boundary condition \-j(0) - 0 . 
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A solution is found by noting that X^ and x^ are equivalent 
in that they have the same differential equations and similar boundary 
conditions: 


(tx^)^+r 2 tx 2 = o 

\(0) = 0 , 


> 

(ex')' h- r 2 tx 1 = o 

x 1 (0) = 0 , 

tx'(l) « T 

« 


Since the differential equations are linear, this requires that 


X 2 = " c lV T ' (5-27) 


where is an unspecified constant. 

Similarly, it can be shown that 


‘i ■ c iV T 


(3.27a) 


Substituting these relations into Eq. (3.26) gives 


2 „2 2 

X 2°l ^ °1 X 1 

1 - ~ = 0 

t T T 


( 3 . 28 ) 


Since -x. = tx^ , this can be written as 


/ \2 


(<> 


P ? 2 

B + r^xf 


(5-29) 


where 


B 2 = T/c x . 
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Murphy (Ref. 58 ) lists three solutions to this differential equa- 
tion, but they are essential Ly equivalent and can be expressed in the 
genera] form 


x^ = ± |r sinh (C ± Ps) . (3*30) 

Here C is the undetermined constant of the differential equa- 
tion. Applying the boundary conditions on x^ gives 

x^ = D sinh Ps/sinh P . (3*31) 

2 p 2 

Note that this determines that T/B = T sinh V/ (Dr) 

Placing this value for x^ in the original differential equation 

gives 

cosh Ps „ sinh Ts dt 

t'FD + 2r t — = - 2r tanh Ps ds . (3.32) 

sir.h T s inh P t 

Integrating both sides: 

in t = - 2 in cosh Ps + , 

or 

t = Cg/cosh 2 Fs . ( 5 * 55 ) 
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The relation tx!^(l) = T provides a value for 
the optimal thickness distribution: 


and hence for 


and 


T 

c p ~ — cosh F sinh F 

TD 


T cosh F sinh V 

t - “ p . (3.3*0 

I'D cosh Fs 


This is the result found by Icerman while including the constraint 
that the first natural frequency must be greater than the excitation 
frequency. This constraint was not explicitly included in the present 
formulation, but it is clear that the constraint is satisfied since the 
solution is identical to Icerman 1 s. 

The question of whether additional solutions exist that do not 
satisfy the frequency constraint, and, if so, what they are, t,ok up 
a large part of the time spent on the thesis. The answer to the ques- 
tion of existence is clearly "yes" and can he demons tra--.ee by looking 
at the behavior of the solution as F becomes large. total weight 

of the structure is proportional to: 

1 

J = J" tds = T sinh 2 r/F 2 D . (3.35) 

0 

J increases monotonically and without limit as T increases. 

The curves of Fig. J.l show that a uniform rod can also satisfy the 
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X (1) = D constraint at any number of excitation frequencies. Clearly 
then,, the uniform rod at some frequency satisfies the constraint and 
has less weight than the "optimal' 1 solution. This indicates that the 
solution of Eqs . (3 -31)~ (3-3*0 is not a global solution for all fre- 
quencies . 

Once this fact is established, the unanswered question is: "What 

are the other optimal solutions?" At first, it was thought that addi- 
tional solutions could be found for Eq. ( 3 . 29 ). After a long fruit- 
less search for other solutions, it was determined that the problem 
was ill-posed, in a special sense. 

The adjective ill-posed has generally been reserved for formula- 
tions that possess no solutions or no physically meaningful ones. A 
structural optimization example of such a problem is that of finding 
the minimum weight thickness distribution for a cantilevered rod with 
the constraint that the first natural frequency of the optimum rod 
have the same naturaL frequency as the uniform rod. If the rod is 
modelled in the same way as was done at the beginning of this sec- 
tion, it is relatively easy to show (Ref. 8) that this problem state- 
ment is satisfied by a uniform rod cf vanishingly small thickness, a 
physically uninteresting solution. 

Since Eq , (3-3*0 gives one solution to the problem at hand, it 
cannot be considered to be ill-posed in a strict sense. However, by 
modelling the rod with three equal length segments, each with constant 
thickness, it is possible to find analytical solutions that satisfy 
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all the boundary conditions and constraints and that have a vanishingly 
small thickness distribution. The diagram below gives a qualitative 
comparison of the mode deflection shape given by Eq. ( 3 . 3 I) and the 
mode shape that this physically unrealistic thickness distribution 
woulu have. 



As the thickness goes to zero for the second solution, the dis- 
placement is unbounded, except for finite values at the root and tip. 

A physically meaningful problem statement must, therefore, have 
additional constraints on the response or involve changes in the sys- 
tem equations themselves. Possible modifications include: 

(1) Imposing a minimum allowable thickness constraint, 

(2) Additions of non-structural mass along the rod. 

(3) More constraints on the response quantities (e.g., inequality 
constraints on the stress or the displacement). 
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The first two modifications were successfully applied to the opti- 
mization problems with natural frequency constraints but have been un- 
successful for the forced case developed above. It is felt that, even 
with a minimum thickness constraint or a non-s tructural mass addition, 
an optimal structure with the frequency of excitation greater than a 
structural natural frequency has discontinuities in thickness. Specif- 
ically, it appears likely that the optimal structures have concentrated 
masses; i.e., thickness distributions that include terms of the form 
t c &(s-s ) , where & is the dirac delta. The motivation for this 

speculation comes from solutions obtained using finite element models 
and piece\\?ise constant continuous models. More comments on this are 
offered at the end of the chapter. 

Inequality constraints, such as those mentioned above, can be in- 
cluded in the manner described in the original formulation. Unfortun- 
ately, the added complexity has made the problems so far insoluble by 
analytical means. As mentioned, there is no reason why the equations 
could not be solved by numerical means. However, once the decision is 
made to go to the computer, the most efficient means of attacking these 
problems is by the use of finite elements. The next sections detail 
how this can be done. 

Before proceeding to this analysis, it should be stressed that 
finding additional function space solutions remains as a suitable goal. 
Variations on the example above are the only analytical solutions for 
harmonically excited structures, as far as is known. Additional analytic 


- ?8 - 


solutions would aid tremendously in uncovering the special features of 
this type of problem. 

D, FINITE ELEMENT SOLUTIONS 

The frustration encountered while dealing with the function space 
formulation led to efforts utilizing finite elements. In any realistic 
problem, the use of finite elements is practically a necessity; but the 
generality and elegance of function space solutions makes them the first 
choice for preliminary investigations. 

Examples are given below that extend the two element case of Sec- 
tion III.B to similar structures modelled by up to ten finite elements. 
Further examples deal with a cantilevered beam structure modelled by 
various numbers of elements. 

The constraints used for these examples are inequality constraints 
on the stress. The Appendix indicates how the stress can be determined 
as a function of the displacements. With this formulation, the aug- 
mented cost function has the form: 

n n 

* - E h - r £ [1 - V*««> Sl • »-5 6 > 

i=l 1=1 

Note that the constraint Js.| ^ S is handled by squaring the stress 
values, thereby obviating the need for absolute value brackets. 
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The thickness is transformed by a technique motivated and de- 
scribed in Section II*C: 


t. 

1 


t . + 

min 


1 2 
2 u j 


( 2 . 11 ) 


The 
the cost 


^ are considered the design variables. Derivatives of 
unction with respect to u. are given by 


a* ( ^ " s aSi/at, 

= u. { 1 + J V — 1 1 — 1 — s- 

du. J (s ) . , [i - (s./s ri 

° j V v max 7 1=1 ' i 7 max 7 


( 3 . 37 ) 


The specific examples given below develop the values for 

SS./^tj . 

1 . Example: Cantilevered Rod 

This section deals with a cantilevered rod excited by a uniformly 
distributed load in torsion. Figure A.l aids in depicting the nature 
of the problem. The steady state equation of motion for the problem 
is given by: 

(- a? [M] + [K])[e} = [p] . (5.38) 

The stresses in the elements are developed in the Appendix: 


GRn 

s _ — q f (Cont’d) 

1 L 1 
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and 


GRn 


[0 i - ®i-i ] 


i — 2, 3? • • < j n 


(A*9) 


By taking the derivative of Eq. (3«3®) with respect to t^ , 
expression for the vector is obtained: 


an 


(- tnf [M] + UK]) 




at j 


aw atK] 


- CO 




C®} * (509) 


*1 


Note that [ 0 } and CcJ^/cJtj} in Eqs . (5*58) and ( 5 . 39) have the 
same coefficient matrix. This fact can be exploited by using a sub- 
routine that solves [A]{x] = {b] by decomposing the [A] matrix 
(Ref. 39)* Since the [A] matrix remains unchanged, it has to be 
decomposed only once to solve for the n + 1 systems of n simul- 
taneous equations 
j = lj 2 , . . n . 

With [^d/^t 
directly: 


1 find 

the separate 

determined, 

the s 

ft , 

GRn 

6 ® 3. 

at j " 

L 

at 3 

ft . 

GRn 

‘ft 

Stj 
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i-1 
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All the tools necessary for a solution using the techniques of 
Chapter II are now assembled. The numerical values used in the com- 
puter program were 


G 

S 

max 

R 

L 


«£ 

3.75 X 10 psi 

J 0 - 

2 ttR 5 * 1352 in 5 

h 

5.3 X 10 psi 

P s “ 

0.1 lbm/in^ 

6 inches 

II 

J§ 

J Q - 4.2 slugs 

120 inches 

P x = 

33,200 in-lbs/in 

t . 
mm 

0.02 inches 



A check on the algorithm was made by first solving the - 0 
case. By using the methods of Section III.C, an exact answer can be 
found for the optimal solution for this statically loaded case. With 
the values of the structural parameters given above, this solution can 
be written as 


p RL 

t = — (1 - s) » 0.3^ (1 - s) . (3.41) 

S max ^0 

Figure 3.3(a) shows a comparison of the optimal solution obtained 
using ten finite elements with the exact analytical solution. The 
agreement is seen to be excellent. 

—2 2 2/2 

Figure 3*5(b) shows a ten element solution for s JW* GJ 0 
“1.0 . It is seen that the effect of the excitation is to make the 
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(b) .T = 1.0 

' ' (2 

FIG. 5*5 — Optimal Thickness Distribution for a Cantilevered 
Rod Using Ten Finite Elements. 
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thickness greater all along the span, compared to the static case. 

This is because the inertial loads act in phase with the excitation, 
necessitating a stronger structure. 

Finally, Fig. ,6 shows two solutions for the case = k.o 
The solution of Fig. 3 .6 {b) is an example where the fundamental frequency 
is less than the excitation frequency and is designated the second sol- 
ution. This second solution is lighter than the first solution by a 
factor of I .36 to 3 * 93 * 

Table "j ,1 compares the rotational displacements and the constraint 
values for these two solutions. The two deflection shapes are seen to 
have similar magnitudes but the second solution is 180 ° out of phase 
from the excitation. This allows the inertial load to partially can- 
cel the effects of the excitation, with the result that much less 
structure is required to satisfy the constraints. These constraints 
are presented in the form g, = [1.0 - (S,/s )^] in Table 3*1* With 

the convergence criterion used for the particular example, a value of 

g. that is less than 0.1 can be considered an active constraint, 'The 
x 

root element of the second solution is at its minimum thickness and 
the constraint is clearly not tight for this element. 

The constraints results for the first solution suggest an inter- 
esting question: "Is the first mode solution fully stressed?" The 

results presented here are ambiguous with the minimum thickness con- 
straint clouding the issue further. One might suppose that it would 
be possible to hypothesize that the optimal solution is fully stressed 
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THICKNESS THICKNESS 





Element 


Displacement 


First Second 

Solution Solution 


Constraint (g^) 


First Second 

Solution Solution 



0.029 


0 . 02 ? 


0.02k 


0.150 


2 

0.058 

- 0.056 

0.012 

0.030 

5 

0.088 

- 0.085 

0.010 

0.016 

4 

0.11? 

- 0.115 

| 

0.006 

0.008 

5 

0.146 

- o.ikk 

0.006 

0.002 

6 

0.176 

- 0.173 

0.006 

o.oxk 

7 

0.205 

- 0.202 

o.ook 

o.okl 

8 

0.23k 

- 0.231 

0.006 

0.018 

9 

0.26k 

- 0.2k2 

0.006 

0.851 

10 

0.293 

- 0.246 

0.006 

0.987 


TABLE 3 , 1- -Properties of the Two Thickness Distributions of Figure 3 . 6 , 

and use the function space methods on Section IIX.C to test the hy- 
pothesis. However;, even for this simple problem,, the analytical com- 
plications make a closed form result impossible. A much simpler means 





















































of testing the hypothesis is available, however. This is the two de- 
sign variable example of Section III.B. Figure 3*3 ( a ) shows an example 
where the first solution is not fully stressed. For this figure, the 
local optimum with the thickness values [t] = {l.27 , O. 35 } has a 
constraint vector given by [0.2h , 0.00} ; i.e., the first element 

is not at the maximum allowable stress while the second is □ 

It should be admitted that the above demonstration is not a rigor- 
ous proof that the optimal continuous structure is not fully stressed 
and that the question merits further study. 

2 . Example: Cantilevered Beam 

The second calculation examines the structural optimization of a 
beam excited transversely by a harmonically oscillating load. As in 
the previous example, a stress constraint is imposed, and it is first 
necessary to derive an expression for the derivative of the stress with 
respect to the design variables. 

The Appendix shows that the stress at the center of the element 
can be expressed as: 


Edn 



Edn 

^ — ( W 21 ** W 2i-2^ i = 2, 37*r»,>n . (A.l8) 
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Physically, this equation says that the stress is proportional 
to the change in the end slopes of the elements. The analysis of 
Eqs. (3.38)-OA>) can be repeated almost directly to give: 


(- of [M] + [K]) 


Sw ) 

1 

2 

S V 

- ' ( 

- or 
e 


Edn 



_ ■ ' 


a h 

2L 

atj 

as i 

Edn 



2 L 

\ St 


5tM] a [K] 
+ ' ' 




at, 


[w] , (3.U2) 




(i — 2, 3 j ♦ • * ) n ) 


The parameters chosen for the optimization program were 


E = 10.5 X 10 psi 

L = length = 120 inches 
d = depth = U inches 
b = width = 12 inches 


t . 
ram 


max 


0.1 lbm/in^ 
100 lbs /in 
0.02 inches 
30,000 psi 


Solutions were found using five elements for excitation frequen- 
cies ranging from *1-2.5 rad/sec to 300 rad/sec. Figure 3.7 shows first 
type of solutions for = h2.5 rad/sec and 1*10 rad/sec . The line 
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THICKNESS (inches) 


"\ 

\ 

X 

X 


WEIGHT = .1754 


EXACT RESULT FOR o> e = 0.0 


.40 .60 

S f-X/L) 

(a) £0 = l| P .5 rad/se.c 



WEIGHT = ,2000 


(k) - 100 rad/sec 

FIG. 3.7-optiMl Thickness Distribution £or a Cantilevered 
Beam Using Five Finite Elements 



superimposed on. Fig, 3*T(a) is the exact solution for the statically 
loaded structure given by Eq. ( 3 . 21). For the parameters given above, 
the exact solution is t -- 0 .5 (1 .0 - 2s + s ) , Even with the har- 

monic excitation, there is close correspondence between the two solu- 
tions , 

Figure 3.8 shows two solutions for t« e =80 rad/sec , The second 
solution is slightly lighter for this case. Another second type of solu- 
tion is shown in Fig. 3«9( a ); while Fig. 3*9(h) plots the weight of the 
two solutions as a function of frequency . It is seen that the first solu- 
tion is the lighter for values of the excitation frequency less than 75 
rad/sec and that the second solution becomes significantly lighter for 
higher excitation values. 

E. CONCLUDING COMMENTS 

Cassis (Ref. 2^) reported on the existence of disjoint feasible 
design spaces in connection with problems dealing with truss structures 
excited by half-wave sine pulses,. It is felt that the problems investi- 
gated in this chapter add a great deal to the understanding of this phe- 
nomenon, primarily because the simplicity of the formulation permits a 
minute examination of the behavior of the structure. The main conclu- 
sion from this investigation is that the natural frequencies play a 
central role in creating the many feasible regions. Structures respond 
vigorously when excited near a natural frequency, accordingly, the op- 
timal designs try to stay away from these resonant conditions. 
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THICKNESS (inches) THICKNESS (inches) 


,40 — 



WEIGHT = 0.190B 



.40 .60 
S 


(a) Structure’s First Natural Frequency Greater Than 


WEIGHT * 0.1633 



(b) Structure's First Natural Frequency Less Than 


2.8 — Two Solutions Cor the Optimal Thickness Distribution 
of a Cantilevered Beam Excited at 'J' ~ 80 rad/sec « 
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FIG. 3.9'a) — Optimal Thickness Distribution for a Cantilevered 

Beam Excited at 1 YO rad/sec (r. > Structure’s First 

Natural Frequency). 



cd e { rad/sec) 


FIG. J ,Ofb'. - -Comparison of the Weights of Local Optima as a Func- 
tion of Excitation Frequency the Structure’s 

First Natural Frequency). 
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The construction of analytical solutions by the methods of Sec- 
tion III.C would further aid in the understanding of these types of 
problems because they show the role that various parameters of the 
problem (such as load, frequency and the constraints) play for a 
range of values rather than the specific values of a particular 
numerical solution. It is currently felt that much of the diffi- 
culty in attaining these analytical solutions is due to the fact that 
they often contain concentrated masses. At the present time, this is 
just a hypothesis that is partially based on the results shown in 
Figs. 5.5(b), 5.8(b) and 5*9( a )* I n these figures, it is seen that 
the elements at the tip are significantly larger than the other ele- 
ments. Based on further studies that used more elements, it appears 
that in the limit as n w the final element is discontinuous from 
the rest of the structure and, in fact, represents a concentrated mass. 
This is an area of current research and efforts to prove (or disprove) 
the hypothesis have so far been unsuccessful. It Is mentioned here to 
indicate the quirks these problems can have and to hopefully aid in 
further research in this area. 
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CHAPTER IV 


WHITE NOISE LOADING 


A. INTRODUCTION 

This chapter moves from the area of the previous chapter, where 
the structure was excited at a single frequency, to cases where the 
structure is excited at all frequencies. In particular, this chapter 
deals with excitations that possess a Gaussian probability density 
function and a power spectrum that has a constant value for all fre- 
quencies. The present analysis considers loads that are random in 
time only. It is possible to conceive of structures that are loaded 
randomly in space as well and of structures whose properties are de- 
scribed in a probabilistic fashion, but these complications are not 
considered here. The motivation for this type of formulation comes 
from the atmospheric turbulence example of the next chapter. The 
turbulence wavelengths are frequently so large that any variation in 
the turbulence magnitude across the span of the wing can be considered 
negligible compared with the time variation due to the aircraft's rapid 
penetration of the gust field. 

The flat power spectrum mentioned above is a useful analytical 
concept and is frequently referred to as a "white noise" spectrum. 

Since the excitation is described in probabilistic terms, it is 


necessary to use probabilistic estimates for the response quantities 
as well. The most useful of these, the mean square values of re- 
sponses, are obtained by integrating the power spectrum of the re- 
sponse over the entire range of frequencies: 

no 

4 - / v°> • c*- 15 

-CO 

It has been shown (Ref. ^0) that the quantities that are of in- 
terest here, the displacements and the stresses, have finite mean 
square values even though the excitation has a finite value over an 
infinite range of frequencies. This fact is very important since it 
allows the development of analyses using the attractively simple white 
noise model. It is, of course, necessary to include structural damp- 
ing in the model in order to obtain a finite response. 

It is not possible to have a disjoint feasible design space for 
this problem. The disjoint properties of the examples in the previous 
chapter arose because of the relationships between the excitation fre- 
quency and the natural frequencies of the structure. Since the white 
noise excites the structure at all frequencies, it is no longer pos- 
sible to have these relationships and, in fact, the design space ap- 
pears to be very well behaved for these problems. The next two sec- 
tions develop the constraint criteria used for the study and the anal- 
ysis needed to evaluate the constraints. These methods are then ap- 
plied to beam and rod models, and optimizations are performed. 
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B. 


FAILURE CRITERIA 


A difficulty intrinsic to the analysis of structures excited by 
random loads Is that explicit values of the response quantities can- 
not be obtained. Instead, mean values or expected values are computed 
using principles from probability. A further complication is that it 
is often unclear what meaning these estimates have relative to the 
safe design of a structure. The aim of this section is to describe 
and evaluate methods that can be used to estimate the life of a struc- 
ture subjected to random loads. 

Cyclic loading, characteristic of white noise excitation, can 
cause a structure to fail even when the magnitude of the applied 
stress is well below the theoretical yield stress of the material 
used. These fatigue failures, which are a common source of failure 
in actual structures, are quite difficult to predict even empirically. 
This is an area of intensive active research that is generally desig- 
nated fracture mechanics. Current efforts divide the fatigue process 
into three separate areas: (l) crack initiation, (2) crack propaga- 

tion, and ( 3 ) strength degradation and failure, A recent summary of 
this type of analysis is given by Yang and Trapp (Ref. 4l), These 
analyses require the definition of parameters relating to load time 
histories, crack size, material properties and other factors, in addi- 
tion to involving lengthy calculations. While the reliability esti- 
mates obtained through the use of these methods should be quite good, 
it is felt that the complexity of the calculations involved makes them 
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ill-suited for the present preliminary analysis. Instead, assumptions 
were made that allowed relatively simple calculations and that re- 
quired the definition of a minimum number of parameters. These as- 
sumptions were obtained from Lin (Ref. 31) with supporting material 
from Powell (Ref. bj) • 

With stochastic excitations, there are two logical failure cri- 
teria, corresponding to tx*o separate modes of failure, that could be 
used in the optimization procedure. The first type is failure due to 
the stress exceeding some specified upper limit. This is commonly 
referred to as first passage or first excursion failure. The other 
type of failure mode treats the damage to the structure as a cumula- 
tive process resulting from the fluctuations in the load. When the 
accumulated damage becomes equal to some specified value, the struc- 
ture is assumed to have failed. (it should be mentioned that while 
this analysis treats these types of failure separately, the more re- 
cent fracture mechanics studies combine these two modes by postulating 
that the random loading causes damage through crack initiation and 
growth which results in the reduction of the failure stress so that 
the final failure is of the first type.) 

The reader's familiarity with certain concepts of probability 
theory is assumed in the following discussion. Papoulis (Ref. ) 
was found to be a useful text for reviewing this theory and should 
aid in the understanding of the pertinent results described below. 
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1. 


First Excursion Failure 




In order to determine an estimate of the time to the arrival of 
the first stress greater than some specified value, it is advantageous 
to make a number of assumptions regarding the nature of the excitation 
process. Basic assumptions are that the process is stationary, Gaussian 
and with a zero mean. If this process is denoted by x(t) , then the 
time derivative of the process, x(t) , is also stationary, Gaussian, 
has a zero mean and is independent of x(t) . The joint probability 
density function and x(t) and x(t) is 


p xA (*,*) 


1 

— * exp 

2ttct cr. 

X X 



(4.2) 


The parameters a and cr. in the above equation are the root mean 
square values of x(t) and x(t) respectively. These can be evalu- 
ated from the power spectrum of # (to) by the formula of Eq. (U . 1 ) : 

XX 


cr 2 = / <J> (cr) dor 

X J XX ' 

~00 

00 

cr? a f cu 2 $ (to) dto 
X J XX 

^00 


P‘- 3 ) 


A second assumption is that large values of x(t) arrive inde- 
pendently of one another. (Ref. shows that this assumption is quite 
conservative for narrow band processes.) This assumption leads to a 
Poisson probability function for the number of times, n , that a 
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large magnitude, U 


Is exceeded in time interval , t 


i 


= — exp (- Xt) . (4.4) 
ni 

The Xt term is the expected number of times the load will be 
exceeded in time interval t . Figure 4.1 helps in explaining this 
and in bringing out a further point. 



FIG. 4.1- -Exceedances of U 

In the diagram, an exceedance occurs when x(t) crosses through 
U with a positive slope or through - U with a negative slope. In- 
cluding the negative exceedances can be justified by the physical argu- 
ment that the examples presented later deal with bending stresses in 
structures that are symmetric about their neutral axis. Therefore, a 
compressive stress of magnitude U is accompanied by a tension stress 
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of magnitude U on the opposite surface. Note that since the process 
has a zero mean, the number of negative exceedances can be assumed 
equal to the number of positive exceedances. 

With this formulation, Eqs . (^.2) and (b.H) provide the basis for 
determining the expected time to the first arrival of value U. The 
X term of Eq. is twice the expected number of positive exceed- 

ances of U per unit time. After placing x(t) = U into Eq. (4.2), 
the expected number of exceedances can be determined by use of the 
formula for expected value: 


X = e(n d ) 


= 2 


/ 


x 


exp 


2tpt cr. 

X X 


(- 1^/20^ - x 2 /2o\) dx 


1 

7T 



(M) 


With the use of Eq. (^.4), the probability of failure in time 
interval t is simply one minus the probability of no failure: 


P F (t) 


1 


-Xt 

e 


(4.6) 


The probability of failure a_t time t is found by differentiating 
Eq. (4.6) with respect to t . The expected time to failure is then 
found by multiplying this probability density function times t and 
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integrating it over times ranging from zero to infinity: 


E(T) * j tX e“ Xt dt . (4.7) 

0 

Integrating by parts yields 

E(T) = l/X . (4.8) 

Equation (4,8) can now be coupled with £n . (4*5) to provide the 
means for determining the constraint on the life of the structure due 
to first excursion failure. If it is specified that the stress in the 
structure cannot exceed some specified value U in the time period 

D 

L , the constraint can be written in the form: 

O 

p p 

gj_ = 1 - Lg'/r — exp (- Uj/2t Tp Z 0 . (4.9) 

°s 

Here cr and cr* are the root mean square values of the stress 
S b 

and the stress rate. 

This constraint is applied independently to each element in the 
structure. It should be mentioned that the concept of fleet or lot 
size has been ignored here. Frequently, first excursion failure is 
defined as the time to failure of just one member of a larger sample. 
If the arrival times of the loads are independent from one sample mem- 
ber to another, the expected time to first failure of one structure in 
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a sample size of n is simply l/nk . This would impose a more 
severe constraint on the individual structure, but, as was mentioned, 
this concept was arbitrarily disregarded. 

2 . Fatigue Failure 

An evaluation of the fatigue life can be made using some of the 
results from the previous section, but it also requires further con- 
cepts. An assumption that makes the fatigue life calculation ana- 
lytically straightforward is one that has come to be known as the 
Pa lmgren -Miner Theory (Ref. U 5 ), This "theory" is based on che 
physically observable fact that a tension specimen that is loaded 
cyclically at a constant amplitude of stress, S , fails in fatigue 

after approximately N cycles. It is postulated that a structure 

b 

that is loaded at this same stress level for q cycles {rr < N ) 

b b b 

has been damaged to the extent that it is at the q Q /N fraction 

b b 

of being failed. It is recognized that experimental results do not 
always support this theory, but it provides a simple general rule 
adaptable to analyses of the type presented here. 

This theory is applied to a continuous random process by deter- 
mining the rate at which peaks of a given magnitude occur. The rate 
of damage is then computed using the formula 

] r q(s) dS 

DR = / , (^S-. 10) 

n N ( s ) 
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where 


ri(S) = number of stress peaks of magnitude 
S occurring per unit time , 

N(S) = number of cycles to failure at stress 
magnitude S 

For the purposes of this work, it is assumed that the damage done 
in a time interval T is simply DR X T . 

The parameter N(s) in Eq . (U.10) can be obtained from curves 
that show the number of cycles to failure as a function of the stress 
amplitude, commonly referred to as S-N diagrams. A convenient ana- 
lytical expression that is used in this work to represent this relation- 
ship, and one that is partially supported by data, is the familiar re- 
lation 


N(S) S b - c . (U-.ll) 

S is the stress amplitude and b and c are positive constants that 
must be determined empirically. This clearly gives N(S) = c/s 

The remaining factor needed for Eq. (U- « XO ) is ri (S ) . Powell 

(Ref, hj) presents an analysis that can be used to readily evaluate 
t\(s) . This analysis starts by modifying Eq. (4.5) to obtain the 

expected number of times a stress exceeds a specified positive value 
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S per unit time: 


E[N + (S)] 


1 

277 



(4.5a) 


The derivative of this expression with respect to S can be 
considered a measure of the number of peaks occurring at the level 
S per unit time: 

aE[N + (S)] S erg 

= 1 exp (- S 2 /2cf) . 

as 2tt cig 

A point that must be considered here is that it is very diffi- 
cult to specify what a cycle is for a random process. Equation (4.12) 
counts only the net number of peaks at level S with the "troughs' 1 of 
magnitude '6 subtracted from the peaks. Figure 4.2 presents the rea- 
soning behind this argument. 



FIG. 4.2 — Peaks in a Record of Random Noise. 


(4.12) 
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There are three peaks in. this diagram at points A , B , and C , 
plus one trough at D * Powell's method says that the damage done by 
this patch of noise is equivalent to the damage done by cycles with the 
magnitude of A , !J , and C minus the damage resulting for a cycle of 
magnitude D . Without belaboring the point, on physical grounds this 
seems to be a better method of counting cycles than one that uses the 
gross number of peaks. Lin (Ref. 42) arrives at the same conclusion 
as that given below by assuming that the process is narrow band. For 
such a process, "troughs" with a positive magnitude are not likely to 
occur so that the problem of net versus gross number of peaks is of no 
importance. Finally, Yang (Ref. 46' derives an expression based on the 
magnitude of the excursion rather than the peak magnitude; this is 
clearly an improvement, but was discovered too late to be included 
in the present work. 

The final step in the derivation is the substitution of the ex- 
pressions for N(s) and q(S) into Eq. (4.10): 

S b+1 or^ exp (- S 2 /20g) dS 
2 cmi^ 

2 / 2 

The integral is evaluated by making the transformation S /2cr s = v , 
leading to 


(fc.13) 



03 

DR = — f (2o- 2 v) b / 2 e” V dv 

n/r n~r ** 


c v 
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This integral can be evaluated by the use of Eq. (3 .$81 ) of 

Ref. 42 


DR 


cr* 


27Tco\. 




i 


where r is the gamma function.. 

To put this in constraint form, it is specified that the struc- 
ture have a fatigue life greater than L f . The constraint is then 
written as 


g g = 1 - DR x L f S 0 . (4.1?) 


This completes the description of the constraints used for the 
randomly loaded structure. It is seen that the structural response 
quantities that are required in order to evaluate the constraints are 
the root mean square values of the stress and the stress rate. The 
next section details how these can be obtained and also develops 
methods for obtaining the derivative quantities that are needed for 
the optimization process. 

G. RESPONSE TO WRITE NOISE 

A finite element representation of the response problem can be 
given by 

[m][w} + [ic]{w] = f{e} . (4.i6) 
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The right-hand side indicates that the equivalent forcing func- 
tion is a scalar multiplying a vector that discretizes the uniform 
load. The scalar F has a white noise power spectrum: 

& FP («0 “ N w - co £ w z ™ .IT) 

Given this representation, the problem is to find the mean square 
values of the stresses, which are in turn a matrix function of the dis- 
placement for the examples dealt with here: 

{S} = [T]{w] . (4.18) 

The exact form of [T] depends on the structure being studied, 
but it is always independent of the excitation frequency and the de- 
sign variables for the present study. 

In order to make the problem meaningful, it is necessary to as- 
sume that the system has damping. Otherwise, the white noise excita- 
tion would result in unbounded resonances and an infinite mean square 
response. This was done by assuming that the structure has damping 
which is manifested by a complex shear modulus or Young's modulus. 

This, in turn, means that the stiffness matrix can be represented by: 

[K] * (l + ia)[K Q I . (4.19) 
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[K q ] is a real matrix that is developed in the Appendix and 
1 + ia is a complex scalar with a representing the damping factor 
which is much less than unjLty. This same representation was used in 
Chapter III and, again. Ref, 36 contains a good discussion of it. 

The response is determined by modal superposition. The modes 
used are the first ran inodes of the system: 

ran 

{wj = 2 a i = > (^.20) 

i=l 

where the [p^]'s are the mode shapes and the a^'s are the modal 
participation factors . The mode shapes are independent of the exci- 
tation while the a^'s are not, so the next step is to determine 
power spectra of the a^'s . 

At a given excitation frequency, , Eq. (4,l6) becomes: 

(- 0 ^ [M] + [K]) [P][a] = F[E] . (4.21) 

T 

By premultiplying Eq. (4.21) by [P] , the equation for [a] can 

be determined as a function of the generalized forces, masses and 
stiffnesses : 


(-co 2 [tn] + PCI) {a} = F[P] T [E] . (4.22) 

3 
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The eigenvectors are normalized so that the generalized masses 
are unity: 


[til] = 

[P] T [M][P] = [ I ] 

> 


IK] - 

[PiVltP] = (1 + id!) [ X ] 


(4.23) 


[ X ] is a diagonal matrix containing the eigenvalues of the system. 

This is a system of ran uncoupled equations that can be solved 
independently for the modal participation factors : 

[- <*>1 + “>?(! + la)] a. - F{ Pl } T {E} a ± 

= fCp.} T Ce}/[- + a£(i + ia)] . (h.2k) 

JL c L 

The term multiplying F is the transfer function H _(jo>) that 

V 

relates to F . This makes it possible to form the power spectra 

for the a's : 


$ 

a. a. 
i J 


H 



(jen) N 


w 


H a 

J 


( 4 . 25 ) 


The bar signifies the complex conjugate. 

The most direct route to attaining the variances of the stress and 
stress rates is to express them in terms of the covariances of the 
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a^’s . For computational purposes, this report distinguishes four 
separate covariance integrals : 


I, = o- 


cl » St 

x i 


/ 


$ do 
a, a, 
i i 


r 2 " Re °a.a. 

1 j 


= Re 


f * 

J Vj 


do 


a. a. 

i i 


/ 


n> <t> ((jo) do 

3< 

1 X 


H = Re Va. 

3- j 


/ 


Re / (d" $ ((d) do 

ct » cl * 

3- J 


(4.26) 


where Re designates that only the reax part is of interest. 

The integrals can be evaluated by making a contour integration 
around the upper half plane. Combining terms from Eqs , (4.17); (4.24) 
and (4.25) into (4.26) gives: 



(4.27) 
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For convenience and clarity^ set (3 - . 

The integrand has no zeroes and four poles : 


z a p(l + ia) 


~ ~ z- 


3 

z k 


= p(i - la ) 5 = , 


- Z-, « - z. 

y 1 


(b.2Q) 


Only poles and z^ are inside the contour. The relation- 
ships between the roots given by Eq. (U.28) and standard contour in- 
tegration give 


I, = C.. 2rrl 


( z i " ^g) ( z x ” s 3^ z i " 


( Z L “ K 1 )( z b " 2 g)( z u " z 5 ) J 


m p 

where = ^w^i * Continuing: 


h ■ a*°i 


2z x 2 fte z^2i <3m 2 fie z^i dta 


7TC, 


2 eftrt z^ | z^ j ' 


(4.29) 
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The magnitude of z ^ is . calculated directly: 


[z y \ S = P 2 (L + o?)= 


The c9ra(z^) calculation, is a bit more difficult; 


z^ = Re(z^) + i c9m{z^) = P (1 + itt) 


By equating the real and imaginary parts of z^ , two 
that can be used to solve for Jm(z^) are formed: 


2 e9m(z^) fte(z^) = pa =*- 

Refz^) = p 2 a /[2 ^1(2^^)] 


[Be(z 1 )] 2 - UM^)] 2 = j3 2 


p\a? 




2 “ [«=5m(» 1 ) 1 


r 1 . 


This results in a quadratic equation in eSm(z^) that 
solution 


,2 , 2 l 


[M^)V 


V 


p ± VP + p a 


(Mo) 


equations 


0 - 31 ) 


. (^-32) 

has the 


(M3) 
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Since c5m(si^) is real,- the minus sign can be rejected and 

JP o 1 JL 

Jta(z-) •= p[(l + a 2 ) 2 - l)] 2 . (4.34) 

1 2 

By substitution of Eqs . (4,34) and (4,30) into the finai result 
of Eq. (4.29) 


C 7T 

I, “ o~E X — o p~T • (4.33) 

1 *J?f3[(l + or) 2 ~ lj 2 p^(l + a 2 ) 2 


Since a < 0.1 ; it is appropriate to make the approximation 


that 


V 


1 + a 


Jl 


a 2 


= 1 + 


+ o(a^) 


(4.36) 


The substitution of the first two terms of Eq. (4.36) into Eq. 
(4.35) gives: 


*1 - 


C^TT 


fpa[i + (a 2 / 2)] 


( 4 , 35 ^) 


It is now possible to neglect the cP/2 term compared with unity 
to get the final result: 


\ = n w (ip.) T £E)) 2 7r/oi?a 


(4.37) 
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The remaining integrals are evaluated in a similar fashion, Since 
the calculations are lengthy, but straightforward, only the final re- 
sults are presented: 


hJp 1 3 t {e3{e} t Cp 3 •net (<J\ + <h. ) 

au i £D j 1 -^i ~ “j ) 2 + C^A)^ +<w j) 2 i 


(H .38) 


S = n w (Cp.} t Ce 3) 2 Tt/^.a 


(^09) 


N w [ P .} T £eHe} t [ p .} 7ja (ta t + (p.) 

2[(p ± - + o).} 2 (a 2 /h)] 


(U-.40) 


The variances of the stresses are obtained by a linear combina- 
tion of the covariances that have just been calculated. The examples 
in the sections to follow use the explicit relationships between the 
stress and the displacement. The general form of Eq. (h.l8) is ade- 
quate for the present derivation: 


[S] = [T]{w} = [T][P]{a] 


The power spectra of the stresses are, therefore, related to the 
power spectra of the modal participation factors by the simple rela- 
tion: 


[* sg (<n)] = [Tns][4 aa (ffl)HP]' E m' 1 . (4A1) 
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The complex conjugate is included in the above equation because 
the [T] matrix contains a complex structural parameter. Since 
neither [T] nor [Pj are functions of the excitation frequency, 
the stress variances are found by replacing $ (o>) with the covar- 

iance matrix for the a's in Eq. (k.ltl), 

tXggl = [T][P][X aa ][P] T [T] T . C'v.te) 

Similarly, 

[X-.] = [T] [P] [X. .] [P] T [T] T 
aS aa 

The square roots of the diagonal elements of [Xg g ] and [Xg*] 
are the rms values of the stresses and stress rates needed in order to 
evaluate Eqs . (^. 9 ) and (4.15). 

It is readily shown that these diagonal elements are real and that 
they involve only the real parts of the [X 3 matrix. To prove this, 
some preliminary notations must be defined. 

Express [T] as (1 + ia) [T Q ] , where [T q ] is real. 

th 

Define tp„ as the i,j element of [Tq][pJ , and pt^ as 
the element of 

Note that tp.. = pt,. and that X „ = X 

ii 11 a. a. a. a 

J 1 3 j i 
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The diagonal elements of the stress covariance matrix can there 
fore be explicitly expressed by: 


mn mn 


v. Sj - (1 - “ £) L ^ E 


a. i 


J-l k-1 

mn r 


lti 


(1 + o?) £ 
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«**«) v a . 
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+ tp. ,tp M X 
JLj ‘ j.j ik a. a, 
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= (1 + 
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a 2 ) S 


L j*l 


(tp. . ) X 
v 1 a, a. 

J J J 


mn 


+ 2 


k=j+l 


tp. .tp. k Re ^ 


(M3) 


All the elements in the equation above are real and, as was to be 
proved, only the real parts of the [X ] matrix are included. This 

1 ao. 

explains why only the real parts of the integrals and 1^ of 

Eq. (^.26) were required. 

This concludes the derivation of the terms needed for the con- 
straint evaluation. A remaining task is the calculation of the de- 
rivatives needed for the gradient in the optimization algorithm. 
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1 . 


Derivative Calculations 


The design variables for these problems, the structural thick- 
nesses, are manifested in the mass and stiffness matrices. The gra- 
dient technique of the optimization algorithm requires that the deriv- 
ative of the constraints be calculated. This in turn requires that the 
derivatives be calculated for all the quantities used to compute the 
constraints and that are a function of the design variables. 

The first step is the calculation of the derivatives of the eigen- 
values and eigenvectors of the system. Fox and Kapoor (Ref. 15) pre- 
sented a straightforward method for calculating these quantities, and 
this method is summarized below. 

Consider the unforced system with a given eigenvalue and eigen- 
vector: 


(- \ [M] + [K] )( Pi 3 = {0} . (4.44) 


For ease of notation, set [F^3 = - X^[M] + [K] . 

The derivative of Eq. (4.44) with respect to the design variable 

t . is 
J 


/ $[M] 

SEK] 
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( 5p i ) 

(- \ 
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- — Em] 

{pJ + ir.] i — } 

l a*. 

atj 




( 4 .^ 5 ) 
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The system given by Eq. (4.44) is self-adjoint so that if Eq. 

T 

(4. 4^) is premultiplied by (p, } , the last term drops out, leaving 



F atKi 

aM 

{?,} [M]{p - } = t P J 

X. 

a _ _ 

at . 
0 j 

L 



[pj • (4.46) 


Since the eigenvectors have been normalized to make the general- 
ized masses equal to unity, the eigenvalue derivative can be expressed 
as : 


* t 

- = {Pj.} 


aw a [m] 

— - v • — [Pjl 

s t j at J 


(4.46a) 


From Eq, (4,44), with the eigenvalue derivative calculated, the 
eigenvector derivative can be solved for; 


— ~ [ Pi ] . 

St, 


(4.47) 


But since [F^] is singular, another equation is needed to specify 

the magnitude of fc)P./dt-3 • This equation comes from differentiating 

1 J 

the generalized mass : 




0 = 2{p i ] J '[M] 


T bW 

+ (p i 3 T C Pi 3 

at. 


(4.48) 
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Equations (4.47) and (4,48) can be combined to give: 


[ V 

N - - 

aIF i ] / at j 

2lp,} T [H] _ 

IstJ 

_ (aiM]/s t.)_ 


CpJ • (4.49) 


In order to obtain a square, non-singular matrix, both sides of 
Eq. (4.49) are premultiplied by [F^ , 2 [M]{p^} ] to obtain. 


EF ± f + 4[M]{p i l[p i } T [M] 


^i 


- [FJ 


2[Ml[p i }(p i } T 




J 


(Pf} 


(4.50) 


This is a matrix equation that can be used to solve for the eigen- 
vector derivatives • Note that since the matrix multiplying 

the eigenvector derivative is not a function of the design variable, it 
is necessary to decompose this matrix only once to solve for the n de- 
sign variable derivatives. A further note is that experience with this 
method has indicated that it is frequently helpful to multiply Eq. (4.48) 
by X as a scaling procedure. 

The remaining steps in the derivative calculation are much less 
complicated. The derivatives of modal covariances and 1^ are 

given below as an example, but it seems of little purpose to show the 
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entire anaLysis here. A few terms must be derived first: 


COT =* X. 
x x 
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1 c)^J 
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This should indicate that the remaining derivative calculations 
are tedious, but uncomplicated. It is mostly a matter of the contin- 
uous application of the chain rule until the final derivatives that 
are required are reached. These are the derivatives of the constraints, 
the first of which is given in Eq. (4.9); 


. t a t -#°§ a n 

gj_ = 1 “ hg — ■ e SO 


The derivative is : 


*1 ■*- 
- = ( Sl - x) — 


i / o: 


°s W 


a S 


■ ( 1 *.5 1 0 


Similarly for , from Eq. (4.15) 


s 2 - 1 - l f -3L ^ r (ill 

2 F 27TCC7 S s \ 2 


And the derivative is 
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D. EXAMPLES 


As in the previous chapter, cantilevered rods and beam examples 
were optimized. Figure 4.3 (a) shows the power spectrum of the white 
noise excitation while Fig. 4.3(b) is a qualitative depiction of a 
response quantity. The peaks on the latter figure represent struc- 
tural resonances which are the main contributors to the mean values 
of the response. 

It is perhaps necessary to justify the use of a finite number of 
modes to represent the response of a structure excited by loads with 
a white noise spectrum. As mentioned in the introduction to this 
chapter, Bogdanoff and Goldberg (Ref. 4o) show that the mean square 
values of the stress and displacement in an Euler-Bernoulli beam are 
finite when the beam is excited by the noise. They do this while 
taking into account an infinite number of modes and by assuming con- 
stant viscous damping. A further indication that a finite number of 
modes suffice is given by Eqs . (4.24) and (4.25) which show that the 
peaks of the spectra for the modal participation factors are inversely 

4 

proportional to . This indicates that the contributions to the 

ruts responses from the separate modes die off quickly as the mode num- 
ber and,, therefore the natural frequency increases. Finally an empir- 
ical justification for using a finite number of modes is given by the 
results below which show that solutions found using four modes differ 
only marginally from solutions using two modes. 
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1 . 


Torsion Rod 


The thin walled rod of Section III.D.l is used again in this sec- 
tion, except that a white noise excitation is now present. The fol- 
lowing list of parameters repeats some of the previous values and adds 
new ones for the special requirements of this problem. 


G = 3-75 X 10 psi 

R = 6 inches 

L = 120 inches 

J Q » 2ttE? = 1352 in 5 

= 12k) (lb) 2 / rad/sec 


p * 0.1 lbm/in^ 

s 

*00 = k - S 3lUg3 
a = 0.05 



U_ ' 40,000 psi 

b 


The parameters b and c are from the equation NS a c and 
were obtained by fitting an S-N curve for aluminum given in Crandall 
and Dahl (Ref. k6, Sec. * The value chosen for Ct is rather 

high and it is recognized that an important part of an actual design 
process using the methods described here would be to obtain more ac- 
curate and justifiable values for the Ct , b and c parameters. 

The constraint placed on the fatigue life was that it be no less 
than one year, and the expected time to stress value Ug was set to 
be no less than one -half year. 


The results of the optimization algorithm are presented in Figs 
U . 1 + and k.5. Figure k.k compares the optimal thickness distributions 
when two, three and eight elements are used to represent the structure 
it is seen that as more elements are used, the total weight remains 
nearly constant while there is some qualitative difference in the dis 
tributions. For the eight element structure, more mass tends to be 
concentrated near the tip. More will be said about this later. 

All the results presented in Fig. k.k used two structural modes 
in their solution. Figure b,$ compares results of analyses using two 
modes and four modes , It is seen that there are some minor differ- 
ences at the tip, but they have to be considered negligible. Table 
k.i gives numerical results for the two cases. 


Element 

Number 

Thickness 



Fatigue Constraint 

2 Modes 

k Modes 

2 Modes 

k Modes 

1 

1 .69k 

I.698 

9.5 * 10' 

5.93 * 10 

2 

1.382 

1 .577 

1.8 * 10 

~k 

3.90 * 10 

3 

1.382 

1.392 

k.o • 10 

l.ok * 10“^ 

k 

1.128 

1.129 

9.1 * 10 ‘ 

2.78 • 10 “ 5 

5 


0 . 873 k 

k .8 • 10" 5 

9.37 • 10 ~ 3 

6 

0.6111 

0.6170 

-4 

5.0 • 10 

l.k 6 * io “ 2 

7 

0.3262 

0.5230 

0.925 

o. 87 k 

8 

- -- 

0.3618 

O.5868 

1.000 

1.000 


TABLE k.l — Comparison of Two and Four Mode Solutions. 


- 105 - 





















THICKNESS (inches) 




This shows that although the four mode solution, took 50$ more 
computer time to converge, it did not appreciably change the results. 

The fatigue constraint values are presented to show that the optimiza- 
tion proceeded to the same level in determining the active constraints. 
The values given are those computed using Eq. (4.15); therefore, the 
numbers near zero indicate that the constraint is almost exactly satis- 
fied (i.e., it is active). 

The optimization seems to have found that placing some weight at 
the tip provides an inertial load that relieves the inboard stress. 

Since this phenomenon is exhibited in the beam results as well, it is 
appropriate to consider this in somewhat more detail, 

2 ,» Effect of a Tip Mass 

This section presents some findings of a brief study that was made 
to justify the optimal solutions that included a large finite thickness 
at the tip. In particular, the study sought to determine what effect a 
concentrated mass at the tip would have on the maximum stress in a can- 
tilevered beam. The hypothesis was that the effect would be to reduce 
the stress. Obviously, this would not be the case for a static loading 
or for a low frequency harmonic excitation, but the results of the op- 
timization indicated that something different was happening for the 
white noise excitation. 

The model studied was a uniform cantilevered beam with a point mass 
at the tip. The excitation was assumed to be uniform across the span 


- 108 - 


and random in time with a white noise power spectra ' density. The mass 
of the beam was kept constant while the tip mass was varied as the only 
independent parameter. 

The problem could be solved by a differential equation approach 
coupled xfith modal superposition as was done in Ref. 4o. However, 
since a computer program that analyzed this type of problem using 
finite elements already existed, it was more expedient to use it. 

The next section presents the structural parameters and the excita- 
tion spectrum used for the analysis. The thickness distribution was 
held fixed for all elements at a value of one. A nonstructural point 
mass was added to the last element and was varied through a range of 
values . 

Figures 4.6 and 4.7 present the results for the rms stress and 
stress rate, respectively, for four values of the concentrated mass, 
nondimensionalized by the mass of the beam. It is seen that the mass 
has the effect of reducing the maximum rms stress, which always occurs 
at the root. The effect on the rms stress rate is to increase its 
peak value, but since the stress is of far more importance in the 
evaluation of fatigue life than the stress rate, this increase is 
relatively unimportant. It is interesting to note that the higher 
modes are obviously present in the stress rate distribution but that 
the first two modes seem to dominate the stress distribution. 

The main finding is that the addition of mass at the tip can im- 
prove the fatigue life. In hindsight it is clear what has happened: 
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RMS STRESS {x 10' d psi) 







the added mass acts as an inertial force that resists the excitation 
and, in the limit as the mass becomes very large, acts as a simply 
supported boundary. 

For the cantilevered rod, a similar effect takes place in that a 
mass would act to restrain the tip rotation and in the limit act as a 
fixed boundary. 

This is an interesting and unanticipated result. A further study 
that could be done is a two design variable optimization study using 
the concentrated mass and the uniform thickness as the variables. 
Constraints could be placed on the rms stress or on the fatigue life. 
The above analysis shows that the optimal concentrated mass would not 
be zero . 

3 . Cantilevered Beam 

A beam example was optimized to see if it had any new, interesting 
characteristics. The methods of Section IV. B are directly applicable 
to the beam example so that the only changes necessary are the inclu- 
sion of the proper forms for the finite element representation of the 
beam structure. Since the Appendix and Chapter III are quite thorough 
in these aspects, they are not repeated here. 

The properties chosen for the beam and the load are 

Length = 2^0 inches E = 10.5 X 10 psi 

Width = JO inches P g = 0.1 lbm/in^ 

Depth = 3*0 inches a = 0.05 
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b 


8 


c 


N w = 0.01 ( lb/ in) /rad/ sec 



10 

4o 


4l 

,000 psi 


The large width to depth ratio was chosen because of a future 
anticipated application of the model to aeroelastic problems where 
it would represent a wing. 

The constraints were continued at one year for the fatigue life 
and one-half year for the ' ■cpected time to failure. 

A comparison of the results obtained using two elements and eight 
elements is presented in Fig. 4.8, while Fig. 4.9 compares results ob- 
tained from an analysis that used four modes with one that used two. 

The concentration of mass near the tip is more pronounced for the prob- 
lem, but the qualitative effects are the same as for the rod example. 

Table 4.2 compares the four mode and the two mode solutions. 


Element 

Thickness 

Fatigue 

Constraint 

First Excursion 
Failure Constraint 

2 Modes 

4 Modes 

2 Modes 

4 Modes 

2 Modes 

4 Modes 

1 

2 

3 

4 

5 

6 

7 

8 

0.1750 

0.1335 

0.101S 

0.0801 

0.0596 

0.0381 

0.0554 

0,644 

0.1751 
0.1352 
0.1033 
0.0809 
0.0602 
0.0387 
0.0590 
0 .0647 

0.042 

0.049 

0.055 

0.069 

0.113 

0.165 

1.00 

1.00 

0.034 

0.0 43 

0.064 

0.070 

O.O96 

0.249 

1.00 

1.00 

0.990 

O.961 

0.941 

0.940 

0.999 

1.00 

1.00 

1.00 

O.99O 

0.977 

O.968 

0.994 

0.999 

1.00 

1.00 

1.00 

Total 

0.7078 

0.7171 



TABLE 4.2--Cantilevered Beam: Comparison of Two and Four Mode Solutions 
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THICKNESS (inches) 



X/L 

FIG, ^.9 — Optimal Thickness Distribution for a Cantilevered Beam Excited by 
a Uniformly Distributed White Noise Transverse Load — Effect of 
the Number of Structural Modes. 


The four mode solution is 1 . 3 $ heavier than, the two mode solution 
a disparity that is probably less than the percentage by which these 
solutions differ from the true optimum. It is possible that further 
iteration would make some of the constraints tighter, but it is felt 
that little Information would be returned to justify the added com- 
puter time. 

E . CONCLUDING COMMENTS 

The results of the two examples tend to show that; as in some of 
the harmonically forced solutions of the previous chapter, there is a 
tendency for some of the mass to be concentrated near the tip. In fact, 
the solutions obtained for the white noise examples could perhaps be 
thought of as a superposition of the two solutions given for a single 
harmonic excitation, such as those of Fig. 3*6* It is not known whether 
this observation has any practical significance for the solution of 
this class of problems . 

It is felt that a formulation of this type makes a useful con- 
tribution in that it presents new results and extends the methods of 
structural optimisation into an almost unexplored field. Obviously, 
however, the examples studied in this chapter are mainly of theoreti- 
cal interest. Methods of fracture mechanics combined with load spec- 
tra that are of great practical interest would aid greatly in the ap- 
plication of the techniques to more applied studies. 
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The next chapter does attempt to perform an optimization on a 
structure that is of more interest: an aircraft wing in the presence 

of atmospheric turbulence. 


CHAPTER V 


CONTINUOUS ATMOSPHERIC TURBULENCE 


A. INTRODUCTION 

Structural fatigue and failure resulting from stochastic loads 
are one of the most commonly occurring maintenance and safety prob- 
lems for aircraft structures. The nature of these vehicles is such 
that there is a very high payoff in terms of performance and operating 
economy for savings made in the structural weight. These two facts 
combine to provide a powerful motivation for finding optimal struc- 
tures under the condition of random aerodynamic excitation with fa- 
tigue life as one of their constraints. Specifically, this chapter 
deals with the minimization of the structural weight of an aircraft 
wing that is subjected to continuous atmospheric turbulence. 

The formulation used in this study is, in keeping with the scope 
of the thesis, of a preliminary nature with a continual tradeoff made 
between physical realism and computational simplicity. The main ob- 
jectives in the development of the mathematical models that are pre- 
sented in the next section are to obtain a representation that is 
consistent in terms of level of sophistication and to retain the 
important elements of the problem. After the presentation of these 
models, it is necessary to develop the analytical tools needed for 
the constraint evaluation and then some results are presented. 
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B. COMPUTATIONAL MODELS 


There are three distinct aieas dint have to be considered in the 
development oil the mathematical representation of a wing excited by 
turbulence: ( 1 ) the structure of the wing, ( 2 ) the aerodynamic oper- 

ators and (j) the disturbing gust forces. Before dealing with each of 
these separately, some general limitations on the analysis should be 
mentioned here. 

The motion of the wing was constrained to consist of rigid body 
plunging motion plus transverse bending, A more general formulation 
would include at least rotational deformation and perhaps rigid body 
rotations as well. While it would not be impossible to include these, 
it is felt that the present formulation is the logical place to start. 

A similar decision was made to limit the constraints to those 
dealing with the life of the structure. It is realized that an ac- 
tual design has to meet a myriad of criteria so that the results pre- 
sented here represent only the specific designs obtained for a specifi- 
cally posed problem. 

1 , Structural Model 

Many of f.he mathematical aspects of the present problem were pro- 
vided by Kef. 32. In selecting a structural model to use in this study 
it seemed natural, therefore, to choose a wing that is used extensively 
in the examples of that text. In particular, Example 10.6 of that text 
presents an analysis that parallels much of what is presented below. 
Figure 5.1 shows a planform of that wing with its important dimensions. 


I I 


^NACELLE 


FUSELAGE 


u 


FIG. 5*1 -“Wing Planform 





As the figure shows, the structural model chosen includes a nacelle 
and fuselage. The masses of these two elements were held fixed during 
the optimization at the values of 

^FUS = ^0.4 slugs , 

ftWc “ l('3*0 slugs 

The assumption regarding linearity between the design variables 
and the structural inertia and mass was retained in this chapter. By 
fitting data given in Ref. 3^, the following factors of proportionality 
were obtained : 


m(y) = mass /inch = 2 t ? t(y) slugs /inch , 

El = stiffness = !5-9^ * 10^° t(y) lbs-ln.^ 

The taper of the chord adds complexity to the numerical calcula- 
tions that determine the mass and stiffness matrices. Section C of 
the Appendix details corrections that are made to the untapered re- 
sults to account for this fact. In addition, the Appendix describes 
how the non-structural masses representing the nacelle and the fuse- 
lage are incorporated into the mass matrix. 


2. Turbulence Model 


The previous chapter dealt with the responses to a random excita- 
tion whose power spectrum was constant over all frequencies. Numerous 
studies have shown that this white noise assumption is inadequate as a 
model for atmospheric turbulence. Chapter 15 of Ref. 47 and Ref. 48 
contain excellent discussions of the procedures used and the approxi- 
mations made in the development of alternative models. From these 
references it was decided that the analytical expression for the tur- 
bulence spectrum that is best suited for the present study is the one 
designated the von Karman model. The power spectrum of the vertical 
component of the atmospheric turbulence given by this model is 

% L T [l + | (1.539 yi) 2 ] 

% (A) = — ^ • (3.4) 

S 7T [ 1 + (1.359 xy) 2 ] 11 / 6 

The terms of this equation are defined in the list of symbols. 

A number of crucial assumptions have to be made about the nature 
of the turbulence in order to arrive at this form (e.g. ? that the tur- 
bulence is homogeneous and that it has a Gaussian distribution) . The 
adequacy of these assumptions are evaluated quite well in Ref. 48 and 
will not be discussed here. 

Values for the turbulence scale and the mean square value of the 
turbulence had to be selected. From Ref. 48, values were chosen that 
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were typical for severe thunderstorm conditions. These were 

l t = 5000 ft. , 

cr =1^ ft, /sec. 
w 

S 

The scale length {which is a measure of the turbulence wavelength) 
is considerably greater than the 85 ft. span of the wing used for this 
study. This large difference in scale reinforces the approximation 
that the turbulence is one-dimensional with a uniform value across 
the span at any instant. 

Figure 5*2 compares the von Kantian spectrum used in the present 
study with the spectrum used in Example 10.6 of Ref. 52 . It is neces- 
sary to present the comparison here because a later figure compares 
two bending moment spectra that were obtained using the two different 
excitation spectra. It is seen that the von Karman spectrum has a 
considerably higher proportion of its energy in the lower frequencies. 

5. Aerodynamic Operators 

The most important difference in the nature of the present problem 
compared to those of the previous chapters is in the manner in which 
the loading is exerted on the structure. In the previous chapter, the 
random disturbance was assumed to be transferred directly to the struc- 
ture in some unspecified manner. In the present example, the aerody- 
namic loads that result from the unsteady gust differ in phase and 
magnitude from the gusts themselves. This is due to the fact that 
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(G)/Vj (ft/rad) 



FIG. 5.2— Comparison of Excitation Spectra. 
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the loads on the structure, which are a function of the circulation, 
do not respond instantaneously to the gust. A further complication 
is the fact that the motion of the wing moving in response to the 
gust's excitation gives rise to additional forces. 

This study is restricted to vertical motions only; therefore, 
the relevant load acting on the wing is the lift. As the previous 
paragraph indicates, this load can be separated into two components: 


P 


L + L 
m g 


(5.2) 


Lg is the direct lift associated with the impingement of the 
gust while is the added lift resulting from the wing's motion. 

Values for these two components are developed in Chap. 5 of Ref. 32 
for a two-dimensional airfoil in incompressible flow that is en- 
countering a sinusoidal guf*- 1 -. These values are given by 


L 


27TP U, 
a d 


c(k)[J o 0O - u-lOO] +iJ 1 (k)[ , 


( 5 - 5 ) 


L = TIP J 2 [k 2 - 2ik C(k)]h 
m a 


(5.>0 


where : 


h 

C(k) 
J 0 and 3 1 


vertical displacement , 

Theodorsen's function , 

Bessel functions of the first kind . 
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Theodorsen' s function is a complex function of the reduced fre- 
quency and is an analytical representation of the change in amplitude 
and phase of the circulatory lift due to a vertical oscillation. It 
can be expressed explicitly in terms of Hankel functions, but for the 
low magnitudes of k of interest to this study, it was deemed ade- 
quate to use an approximation that is given by Fung in Section 6.9 of 
Ref. 49: 

O.I 65 O .535 

C(k) = 1.0 . ( 5 . 5 ) 

[ 1.0 - (0.0 J i55i/k)] [ 1.0 - (0.31 A )] 

Perhaps it is in order to point out here that the complex nature 
of the aerodynamics makes it unnecessary to include structural damp- 
ing. This damping was required in the previous chapter in order to 
obtain finite response, but the out of phase component of the aero- 
dynamics acts as a damping mechanism that limits the structural re- 
sponse to finite values regardless of the excitation frequency. 

In order to apply these results to the problem at hand, a number 
of additional assumptions must be made. These are mainly the approx- 
imations that are used in aerodynamic strip theory: 

(1) The incompressible results are valid for the analysis. (The 
example considered has a free stream Mach number of 0.62 so that com- 
pressibility effects could be constructively considered.) 

(2) The reduced frequency is computed using a reference chord, 
as opposed to the local chord, resulting in a k that is constant 
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across the span. While this is not strictly necessary, it greatly 
simplifies the calculation. The range of actual reduced frequency 
values across the span is small enough so that the error introduced 
by this assumption is not large. 

( 3 ) The loads an the three-dimensional wing arc the same as 
would occur at that wing station in a two-dimensional flow (except 
for the disparity in k values mentioned in the previous assumption). 

It would be interesting, and not too difficult, to determine what 
effect these assumptions have on the final results. However, these 
were considered to be secondary matters that did not require evalu- 
ation for the present study. 

Once these aerodynamic loads have been evaluated, it is necessary 
to put them into a form consistent with the finite element models de- 
veloped for the mass and stiffness matrices. Again, the Appendix pro- 
vides the details of how this is done. 

Finally, values of the parameters necessary for calculating the 
aerodynamic loads are: 

U = 696.8 ft/sec f 

p = 2.378 x 10 ^ slugs/ft^ , 

ci 

b _ = 6.771 ft 

ref 

* - 
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C, RESPONSE QUANTITIES AND GRADIENT EVALUATION 

The end result of the development of the models in the previous 
section is the construction of an equation of motion in the form: 

(- [M] + [K] - [A]) {w} = [G] . (5-6) 

Some new terms have been added to the formulation used to study 
white noise. These are 

{g} = Vector representing the load due to a unit sinu- 

soidal gust of frequency 0^ 

[A] = Matrix relating the loads on the aircraft due to 

the aircraft's oscillation at frequency 

'Hie general method used in the previous chapter can be repeated 
here to find the root mean square response values for the stresses and 
the stress races. However, the new elements of the problem necessi- 
tate going through a brief description of these methods. While it is 
not explicitly emphasized, it must be remembered that the analysis 
presented below is in terms of a unit gust excitation. 

Modal superposition can again be used to obtain the response of 
the wing at a specified frequency: * 

mn 

[w] = 2 a i = ^ ^ * (5-7) 

i=l 
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the [p^} vectors are the eigenvectors of the system 

(- X,[M] + [K]){p,] = 0.0 and the a.'s are the modal participation 

factors that are to be determined for the forced response. The next 

T 

step is to premultiply Eq. (5*8) b y IP] : 

(- a? [ I 1 + [ X ] - [GA]) [a] = {GG} . (5.8) 

The mode shapes have been normalized so that the generalized 
masses are unity. The new terms of Eq. (5*8) are clearly 

[GA] = [P] T [A][P] 

{GG} = [p] t [g] 

In the previous chapter, multiplying the equation of motion by 
the transposed eigenvector matrix uncoupled the equations in the a^'s 
by diagonalizing the mass and stiffness matrices. The generalized 
aerodynamic matrix is not diagonal, however, so the system of equa- 
tions for the {a} vector have to be treated simultaneously. 

Also, since [GA] and [gg] vary in a complex fashion with the 
reduced frequency, it is necessary to evaluate Eq. (5*8) at a number 
of discrete reduced frequency values. 

Once the modal participation factors have been found for a large 
enough number of reduced frequency values to represent the complete 
range of interest, it is possible to move on to the calculation of 
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the stresses. Once again,, the methods of the previous chapter are 
inadequate for this problem. Hie difficulty now is that since the 
model permits rigid body motions, the bending moments cannot be cal- 
culated from a derivative of the displacement vector. Instead, ex- 
ternal and inertial loads are summed and the bending moment is found 
from these forces and from the fact that the shear force and bending 
moment at the wing’s tip are zero. The force acting is given by 

2 

F = L + L + mOi w 
g m e 

Or in matrix notation: 

{F} = ([A] + o£ [M]) jw} + {G} . (5*9) 

The [f] vector represents concentrated forces and moments acting 
at the node points. From this vector, it is possible to calculate the 
bending moment acting at any specified location on the wing. For the 
purposes of this analysis, the bending moments were computed at the 
center of each element plus an additional calculation at the wing's 
root . 

Performing the moment summations at these points gives: 

n r 1 

BM t ■ Z P 2j + 1 + k P 2j + 1] ( C ° nt ’ d > 
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and 


BM Koot = 



F, 


,1. p 

y-ti 2j 



( 5 . 10 ) 


This can be summarized by a matrix equation: [EM] - [T][f] 

The vector of bending moments calculated in this way can be 
thought of as the admittance functions for the structure. The fac- 
tor that is of prime interest is the mean square bending stress. 

Given the bending moment, the remainder of the calculation is quite 
s traight f orwar d. First, the admittance of the bending stress is 
calculated using the standard S - Mc/l formula. Proper account 
has to be made of the tapered property of the wing in this calcula- 
tion as it enters into both the c and I terms in the stress equa- 
tion. With the bending stress admittance calculated at a number of 
frequencies, the mean square response is calculated from: 
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w x ' 
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( 5 . 11 ) 


The mean square stress rate is computed in a similar fashion: 

CO 

4* “ f 0)2 I s I 2 * <P) * (5-12) 

0 g 

Simpson's rule was used in performing the numerical integrations. 
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Once these two parameters have been determined, the analysis of 
Section IV. B can be used to determine the fatigue life and time to 
first excursion failure at the wing stations of interest. This anal- 
ysis will not be repeated here. 

The changes in formulation described above also create some dif- 
ferences in the way the gradients are calculated. Again, only the new 
details are described in this section, since the previous chapter is 
available to provide added detail. 

The eigenvalue and eigenvector derivatives are found in the same 

manner as previously except that the rigid body mode allows certain 

simplifications. Specifically, since the rigid body frequency is zero, 

the derivatives of this frequency are trivally zero: 

1 

d X l 

— - =0 j = 1,2, ...,n . (5*1?) 



The rigid body mode shape is a vector given by relation: 

[ Pi } T = q(l, 1,0,1,. ..,1,0} = ri{u} T , (5-14) 

X 

where q is the normalizing factor used to obtain {p^} = 

1.0 = q 2 {u} T [M]{u} . 

Since the mass matrix varies with the design variables, the 
rigid body mode does have a derivative with respect to the design 
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variables that can be evaluated by the use of the relationship just 
obtained: 


<p Bn 

2n[u} T [M]{u3 — • 

At, 
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' T - — [u] 


sc : 


Bn , T slMl 

— - - n 5 {u] T {u}/2 


(5.15) 
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The matrix triple product [u] (BM/Btj) [U] can be shown to 

til 

be equal to the structural mass of the j element, m. , divided 
by the design variable t^ . The derivative expression for the mode 
shape then becomes : 
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The next step is the determination of the derivative of the modal 
participation factors. Recall Eq . (5*8): 


(• ^ [I] + [ M " [GAJ) {a} = [GG] 


Taking the derivative with respect to the thickness of the j 
element gives : 
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(5.16) 


/a r atn \ , , 

[ [ X ] - rll'PW {a} 

\a*j a c j / 

As in previous cases of this ty^e, the matrices 'on the, left-hand 
sides of Eqs . ( 5 * 8 ) and (5.16) are -he same, regardless of which de- 
sign vector is of interest. Therefore, the matrix decomposition of 
o 

(-0> [I]+[X]- [GA]) needs U) he evaluated only once for the 

n + 1 systems of 2n + 1 simul t artivms equations. 

Another note is that the derivative of the generalized aerody- 
namics matrix involves only the mode shapes since the aerodynamics 
matrix, [A] , is not a function of the design variable. This is 

different from flutter optimization problems, where the aerodynamics 
are indirectly a function of the des Lgn variable because the flutter 
frequency is contained in the matrix (Kef. 5 ^)* 

Finally, note that even though matrix [GA] is symmetric, the 
derivative ^[GA]/^t^ is not. 

The remaining derivative, culcui itlons can now be evaluated: 
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Another new derivative that must be evaluated is: 

S J 2 - as i 

“ — =■• 2S — , (5.18) 

J 1 ° J 

where the bar indicates the complex conjugate. 

Since the bending stress is proportional to the bending moment 
and inversely proportional to the element thickness, (S^ = cp^BM^/t^ , 
where cp^ is the constant of proportionality), the bending stress 
derivative is given by 
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where o. . is the Kronecker delta, 
ij 

Finally, 
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(5.19) 


( 5 . 20 ) 


and similarly for the stress rate. The remaining derivatives for the 
constraints and the objective function are identical in form to those 
of the previous chapter and are not repeated here. 


D. RESULTS 

As the above descriptions have perhaps indicated, the function 
evaluation and gradient calculation require a considerable amount of 
computation. Consider an example that has N elements and a mesh of 
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NF discrete frequencies used in. the response calculations. Further 
specify that FIN natural modes are used for modal superposition. 

Then each function evaluation requires the solution of a 2N + 1 
eigenvalue problem. In addition, the MN linear simultaneous 
equations given by Eq. (5.3) must be solved NF separate times. 

If gradient information is desired, the MN simultaneous equations 
given by Eq. (5*15) must be solved NF x N times. An additional 
factor is that unless one is very clever or sacrifices programming 
speed and clarity, the arrays needed for the computation quickly fill 
the comput c's available core, E.g., a reasonable way to dimension 
^s/^tj of Eq. (5.19) i- s DS(N + 1,NF, N) signifying that each of the 
N + 1 stress values for each of the NF frequencies has derivatives 
with respect to N different design variables. 

For these reasons, the examples done for the thesis were kept as 
simple as possible while retaining the capability of obtaining mean- 
ingful results. 

The first example used three structural elements and retained the 
rigid body mode plus one bending mode. Twenty-nine reduced frequency 
values ranging, at equal intervals, from 0.0 to 0,28 were used. Al- 
though this first example was worked mainly as a check on the algorithm, 
the results are of sufficient interest to be presented here. The con- 
straints were identical with those of the previous chapter in that the 
fatigue life was specified to be greater than one year while the time 
to first excursion failure was specified to be greater than one-half 
year. 
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are 


* 

The initial and optimal thickness distributions for this example 



' 1.00 
* o . 90 ’ 
. 0.50 J 


i 0.04965 ) 
0.02539 j 
0.01146 / 


A plot of the final thickness distribution is given in Fig. 5.3. 
The active constraints designated on the figure are all first excur- 
sion failure type constraints. 

The marked reduction in weight is partially due to the fact that 
the initial configuration is extremely overdesigned with respect to 
the constraints considered here. The rms stress at the root for the 
initial design is approximately 600 psi; a value so far below the 
specified ultimate strength level of 40000 psi as to be insignificant. 
TIi is should not be too surprising, since the textbook example from 
which the model was obtained was not intended to be near a critical 
value with respect to this particular constraint. It is surprising 
that the weight is reduced by a factor greater than twenty. This fact 
is discussed following the presentation of the second and last example. 

The final example used five elements to represent the structure 
and retained two bending modes plus the rigid body mode. The same 


* 

Hie cost function used for the wing examples was the sum of the 
design variables. Due to the taper of the wing, this is not exactly 
proportional to the structural weight; therefore, the final thickness 
distribution is not the minimum weight solution. This oversight was 
detected after the two examples were completed, and it did not seem 
a large enough error to require rc-optimizations with their attendant 
computer costs. 
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THICKNESS (inches) 




ciumber o£ reduced frequencies were used. Since the previous example 
had permitted an. optimal design that was unrealistically light, the 
constraint lives were multiplied by a factor of ten. The fatigue life 
was therefore constrained to be greater than ten years and the time to 
first excursion failure was required to be greater than five years. 

The remaining parameters x?ere left unchanged. 

Figure 5.H compares the power spectral density of the root bending 
moment obtained from Example 10.6 of Ref. J2 with that obtained using 
the initial design and the models developed for the present study. 

The different turbulence spectra used for the two cases account for 
the majority of the discrepancy, while some differences in the model- 
ling of the structure account for the shift in the location of the 
second peak. In the figure, the first peak is almost entirely due to 
the rigid Dody response while the second peak occurs very close to 
the natural frequency of the first bending mode. The second bending 
mode occurs at such a high frequency that it does not have an effect 
on the root bending moment. Of interest here is the fact that the 
two solutions are qualitatively the same, indicating that the computer 
analysis has been done correctly. 

The initial and optimal thickness values, as well as the rms 
values of the stress and stress rate for the final design are given 
in Table 5*1* 

Figure 5.5 plots the final thickness distribution and indicate;; 
where the constraints are active. In this example, all the active 
constraints were of the fatigue type. 
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Element 

Thickness 

Final Design 

Initial 

Optimal 

RMS Stress 

RMS Stress Rate 

Root 

— 

! 

4j558 psi 

25,706 psi/sec 

1 

1.00” 

0 . 0628 ” 

5965 

18,562 

2 

0.955 

0.0312 

5521 

31,506 

3 

0.915 

0 . 02 U 5 

5733 

24,971 

4 

0.810 

0.0125 

5594 

34,877 

5 

0 .380 

0.0035 

5157 

32,590 


TABLE 5.1 — Wing in a Turbulent Atmosphere 

The most striking fact that this solution exhibits is that, even 
with the constraint lives multiplied by ten from the previous example, 
there is a very large weight decrease from the initial to the final de- 
sign, Part, of the explanation for this behavior is indicated by Fig. 5»^. 
This figure shows the power spectral density of the root bending moment 
for the final design. A comparison of this figure with Fig. 5-^ points 
oul two things: (1) The area under the power spectrum, and hence the 
ms bending moment, for the final design is substantially less than the 
area under the comparable curve for the original design and ( 2 ) the re- 
sponse to the first bending mode has disappeared in the optimal design. 
These two results are related to the fact that, as the weight is re- 
duced, the inertial loads become increasingly less important compared 
to the aerodynamic loads . 
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Of course the rms stresses are considerably greater for the final 
design since the stresses are inversely proportional to the design var 
iables . 

This example dramatically illustrates a tenet of structural op- 
timization that is frequently ignored; viz., for the final result to 
be useful, all the design conditions that the structure will be re- 
quired to meet must be considered simultaneously. Because this work 
was primarily interested in studying the effect of stochastic loads 
in the optimization process, other constraints that would have made 
the final design more meaningful were left out of the analysis. The 
following, concluding chapter offers some suggestions as to how the 
optimal design of a wing could include more complete design conditions 



CHAPTER VI 


CONCLUSIONS AND RECOMMENDATIONS 

A number of comments about the behavior and significance of the 
solutions have already been made. Sections II. D, III.E and IV, E are 
devoted to discussions of the material presented in their respective 
chapters. This final chapter reiterates and expands on these comments 
and lists areas that would benefit from further study. 

A general conclusion is that despite the complications introduced 
by dynamic loadings , methods of mathematical programming can be applied 
to studies of this type. This is not a surprising, or even new, con- 
clusion. Accordingly, the main contributions of this thesis must re- 
side in the formulation of the analyses and in some of the interesting 
results obtained. 

In particular, the discovery o f disjoint feasible regions in Chap- 
ter III is of theoretical and perhaps practical interest. The two de- 
sign variable example of Section III. A demonstrates the concepts of 
disjointness with a simplicity that makes it valuable as an instruc- 
tive tool. The ease of the formulation, coupled with the large amount 
of information garnered, also reinforces the maxim that simple cases 
should be examined first. It is felt that some of the studies of op- 
timization with dynamic loading have ignored this rule and have thereby 
suffered from a lack of understanding of the basic principles involved. 


It is realized that the optimal structures shown in Chapter III 
with a first natural frequency that is ]ess than the excitation fre- 
quency are impractical because of their inability to sustain even, 
moderate static loads. However, when a design is influenced by a 
harmonic loading, it may be possible to obtain a more suitable struc- 
ture by "loosening 1 ' it so that one or more natural frequencies are 
less than the forcing frequency rather than by stiffening it so that 
all natural frequencies are higher. 

A design procedure that is related to this concept is that of 
"detuning". In this procedure, masses are added or moved on the 
structure in a manner that minimizes resonant responses . Another 
example is the "soft mounting" of nacelles. In the latter technique, 
the mounts minimize nacelle motion by giving the nacelle-mount com- 
bination a first natural frequency that is less than the predominant 
excitation frequency. Methods developed in this thesis can aid these 
design practices by performing them in a more systematic fashion. 

The results of Chapter IV indicate that the optimal structure can 
differ markedly from intuitive designs that are based on static strength 
alone. Another conclusion is that the response to white noise can be 
adequately estimated by a very small number of modes. For the prob- 
lems studied, the first two modes were of primary importance and further 
modes added little to the final results while increasing the computer 
time for solution significantly. 

The results of the preceding chapter, dealing with a wing in atmos- 
pheric turbulence, make it clear that care must be taken to adequately 
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formulate problems of this nature, The example chosen had an initial 
design that was much too stiff to make a comparison between the ini- 
tial and final designs meaningful. Also a number of features should 
be added to the model studied. Among these are additional loading and 
constraint conditions to insure that the optimal structure has adequate 
strength to withstand normal flight loads and landing impacts. The 
torsional modes are probably also of importance and should be included 
in further investigations. When other factors, such as fracture tough- 
ness, flaw growth, realistic aircraft structures and improved aerody- 
namics are added to this list, the analysis is clearly one that is 
beyond the scope of this developmental work. Hopefully, it has pro- 
vided some of the techniques that a more ambitious research group could 
build upon. 

One aspect of most papers on structural optimization which is ab- 
sent from the present study is the presentation of the amount of com- 
puter time necessary for a problem's solution. Since the time to sol- 
ution was not considered to be an important factor in the problems worked 
here, no attempt was made to minimize it. An indication of the magnitude 
of the computation time required is given by the fact that each design 
iteration for the problem of Chapter V required approximately a minute 
of CPU time on Stanford's IBM 360/67 in the "Quick" partition. By using 
a more efficient compiler and by using approximate techniques for the 
analysis, it seems quite conceivable that this figure could be reduced 
by roughly a factor of ten. For problems with a larger number of design 
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variables, these efficiencies are clearly needed in order to make the 
optimizations manageable from a computer resources standpoint. 

Aside from the above comments on problem formulation and compu- 
tational efficiency, further investigations could be conducted in a 
number of areas. Some of these that the author finds intriguing and 
of importance to gaining an understanding of the principles of struc- 
tural optimization with dynamic loading include the following: 

(1) Further work on the function space solutions of Chapter HI. 
The most desirable goal would be analytical solutions for a range of 
constraint and loading conditions. Unfortunately, analytical results 
to date have been minimal despite considerable efforts made in a search 
for such results. 

(2) Solving the two point boundary value problems of Section III.C 
by numerical methods could also be a worthwhile activity. Pierson (lief. 
10) describes an excellent algorithm that can be used to iteratively 
solve these types of problems. 

(3) Design of optimal two-dimensional structures, such as plates 
and shells, with a harmonic or other dynamic excitation. To the author's 
knowledge, this is uncharted territory and should provide a wealth of 
problems and new results. It is anticipated that the disjoint feasible 
regions will continue to be a complicating factor in the search for op- 
timal solutions for the harmonically loaded structures. 

(^) For problems with stochastic excitation, this work uses rela- 
tively simple failure criteria. The discipline of fracture mechanics 
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has recently developed more sophisticated methods for predicting a 
structure's damage due to random loads. These methods should have 
application to the optimization problems, 

(5) The problems of Chapter IV should be solved for a range of 
such parameters as structural damping, load magnitude, and constraints . 
This would show the qualitative effects that these parameters have on 
the solutions and might even point out some new properties of s tochas tic 
optimization problems. 

These are suggestions that relate directly to the investigations 
of this thesis. There are, of course, many other problems involving 
optimization of structures under dynamic excitation that are of interest 
and importance. It is the author's view that the most critical current 
task for the optimizer is that of acquainting the designer with the tech- 
niques of optimization so that they can jointly determine where the methods 
are of most value. It is felt that this work has broadened the applica- 
tion of optimization methods and has therefore enhanced their attractive- 
ness and usefulness . Hopefully, this promise will be furthered and ful- 
filled . 
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APPENDIX: FINITE ELEMENT MODELS 


The finite element formulation. 1 ; used in the thes is are presented 
in this appendix. The first two sections deal with the representa- 
tion of rod and beam elements with constant cross sections. A stan- 
dard text on finite elements (e.g., Ref. |?l) contains most of the re- 
sults shown in these sections; they are included here for completeness 
and to demonstrate the notation used. A third section deals with the 
representation of the tapered wing used in the fifth chapter of the 
thesis. It is necessary to go into added detail in order to indicate 
the adjustments that the tapered elements require. Lastly, the aero- 
dynamic matrices needed in CahpLcr V arc formulated. 

Throughout the Appendix, a distinction is made between matrices 
that represent a single element and those that represent an assembled 
structure. The former are denoted by a subscript that defines the ele- 
ment being considered (e.g., [K]^ } while the latter have no subscript 

(e.g., tK] ). 

A. TORSION ROD 

Figure A.l shows the rod and the degrees of freedom used. For a 
single element, the mass and stiffness matrices are represented by: 
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( a . Representation Using n Elements 


ORIGINAL PAGE IS 
OF POOR QUALITY 



FIG, A. 1 — Cantilevered Rod. 
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From these elements matrices, the final matrices are readily as- 
sembled. Some assumptions made when this .is done are: 

(1) The one-dimensional structure is divided into n elements 

of equal length: £ - L/n 

(2) The rod is a thin walled tube with structural properties that 
are proportional to the thickness: 


I'GJ) . 

• i 


i'I ). 

1 cri 


GJ o fc i 




(A.2J 


GJq and 1 are constants for the structure. 

(3) The cantilevered boundary condition is accounted for by 
deleting the degree of freedom associated with the root 
station. 
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The matrix is then assembled by adding the contributions of the 
individual elements at the node points. The assembled [M] and [K ] 
matrices have dimensions n X n . 

The equivalent force vectors are needed to represent a uni form 
load across the rod'. They are found by using the formula 

1 

[P } 5:5 p x £ f ^ ds * ( A -3) 

1 0 

Here [a] gives the relation between the discrete displacement 

representation and the continuous actual displacement. For rod ele- 
T 

ments, [a] is given by 

[a} P [l.O ** s , s] . (A.^) 


Then 


{p } 

eq , 


p i 
x 


l. 

2 

,1 

2 


The assembled equivalent force vector is 


{P} = 


P l 
x 


(A. 5) 


(a.6) 
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Another finite clement matrix needed is one for determining the 
stresses. The relation between the stress and the displacement is 
found using: 


where 


s, . o{b] T (e] 


(A.7) 


T s d(a) T R 

fb] T = - [-1,1] 

l ds l 


(A.3) 


Therefore, 


and 


GR 


GR 

S 1 - 7 < e i - 


(A.9) 


B. CANTILEVERED BEAM IN BENDING 

Figure A. 2 shows the beam and the degrees of freedom used in this 
thesis. Note that a slope as well as a vertical deflection are used 
at each node. The mass and stiffness matrices for a single element 
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are given by (Ref. ^1) 




s 




where 


(A. 10) 


[KG] 


[MO] 


~ 12 6i 

-12 

61 

4i 2 

-6i 

2 £ 


12 

- 6* 

Symmetric 


) 1 n 

1 — 


1 1 0 

_ 156 22 l 

54 

-1 3i 

4f 2 

13 i 

- 3i 
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-22 £ 

Symmetric 


4i 


(A. 11) 


As in, the previous section, the structure is divided into n 


equal elements of length t ~ l/n . The inertia and cross sectional 


areas are expressed as linear functions of the element thickness: 


(El). ^ El t . , 

' i Ox J 

(PA). « pA 0 t. . (A. 12) 


The cantilevered boundary condition requires that the degrees of 
freedom corresponding to the root displacement and slope be eliminated. 
Tlie assembled mass and stiffness matrices can again be formed by adding 
contributions from the individual elements at the nodes. The assembled 
matrices have dimensions 2n x 2n . 

The equivalent force vector and the stress vector are found in the 
same way as in the previous section with the {aj and [b] vectors 
replaced by 


{a}' 


{b} 


d 

2.f 


2 

3s + 2s* , 

(s - 2s 2 + .3)1 f 


2 5 

3s - 2s p , 

C-3 2 + S 3 )i3 T f 

(A. 13) 

/ .2 \ T 
(da 

Us S ) 

} 

(A. 14) 

the beam. 



equivalent 

force vector is 


| 

1 \ 


\ 

=* P x * < 
i x 

j -e/12 / 

> 1 ’ 

(A. 15) 

( 

-f/l2 | 
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The assembled vector has the form: 



The stress vector is found using: 


(A,l6) 


=: E[b} T {w] 


Ed 

— — {-6+12s ) (-^H-Gs),* } 6-12s 

2 1 


(-12+6s)f} T 


1 w 2i -3 
w 2i -2 
W 2i -1 
w 2i 


(A. 17) 


If the stresses are calculated at s = ^ , the above relation is 

r T 

greatly simplified in that the first and third elements of {b} are 
zero. This x^as done for this thesis giving an assembled stress vector 
of the form 


(Si 


Ed 
2 1 


0 10 00 - - - 
0-10 100 - - 
0 00 -1 01 - - 


0 00 - ---1 0 
0 00 - - - 00 



(A. 18) 


- 165 - 


C. TAPERED WING 


This section develops the matrices needed for Eq. (5.6): 

(- a? [M] + [K] - [A]) [wj = {G} . (5.6) 

The wing used for this analysis is shown in Fig. 5-1 • The linear 
taper of the chord adds complexity to the analyses of the previous sec- 
tions and these complications are outlined in this section. 

The first step is to represent the nondimensional semi-chord length 
as 


a(y) - b/b 


ref 


root 


ref 


(1.0 - y/tj 


(A. 19) 


For the wing used, the numerical val-ues of the parameters are 
given by 


root 


112.5" j t f - 900" , 


ref 


= 81.2V' 


It is assumed that the thickness to chord ratio remains constant 
across the span so that the aerodynamic thickness is also linearly 
tapered. In order to avoid further complications of doubtful utility, 
the structural thickness is held constant across the length of an ele- 
ment. ’The consequences of these two assumptions are that the mass 
ries as (1.0 - y/t f ) while the inertia varies as (1.0 - y/t-) 1 ^ , 


va 
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As In the previous section, submatrices that deal with a single 
element are the building blorks that make up the final assembled ma- 


trices. It is therefore necessary to have structural properties ex- 
pressed in terms of each element. As an example, the dimensional chord 
til 

length for the i element can be expressed as : 

^(s) - G root Q i (1 " s / 5 i) O^s^l . (A.20) 

By the use of Eq, (A. 19); it can be. shown that 

(X = 1 - L(i - l)/nt f , 

nt nt 

5. = — -+1-1 = — - a. . (A .21) 

1 L L X 


1 , Mass and Stiffness Matrices 

The mass matrix for an element Is found by evaluating the integral 


[M]. 

fp 

where [q] 
tion gives 


= Vs 


A root t f a i f 1 ' t) ™ tnl ' 


ds 


(A. 22) 


X 

is identical to {a} given in Eq. (A.l). The integra- 



te A LOJ. 
x s root i 


n 



(A. 23) 
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[MO] is given in Eq . (A. 11) and 


[Ml] 


72 1 hi 54 -12£ 

3sP iht ~ 5Z 2 

2k> -30 1 

Symmetric p 

5r 


Similarly,, the stiffness matrix is evaluated using 

[Kl i - E {if f h w !p}1 ds 

0 


(A, 24) 


(a.25) 


where 




T 


and 


I. 

l 


I 0 t i l 


.o? 



Since s/&. « 1 , it is reasonable to make the approximation 


1 - 


&. 

i 


l - 2S. 

5. 

l 
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Then, Eq. (A. 25) is evaluated to give 


Ik] 



(a. 26 ) 


[KO] is given in Eq. (A. 11) and 


[Kl] 


6 2 l -6 kt 

2 ? 

r - 2 1 r 

6 -hi 

Symmetric 0 

3^ „ 


(A. 27) 


Note that as t^ , and therefore 5^ , become large, the mass 

and stiffness matrices approach the untapered results of the previous 
section. 

2 . Aerodynamic Matrices 

The aerodynamic matrix [A] results from the motion and can be 
formulated in a manner very similar to that used for the mass matrix. 
From Eq, (5 A), the distributed lift resulting from the motion is given 
by: 


L (v) 
m '* ' 


7TP U 
a 


[a(y)tc] - 2ia(y)C(k) 


(5.ka) 
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Thu equivalent matrix for the iinite element representation is 


found using 


[A]. = 7TP 


f t f 


(a.k)‘ 


i - — 

a. p. 


2ikC (k) Gh I 1 - £-* 


f T ds 


(A. 28) 


p 

Using the approximation that (1.0 -s/h £ )~ : 1.0-2s/fc : 


[A], 


.,2 

7TP U £ 
a 

420 


ak 2 - 2ia kC(k) 


[MO] 


2(a?k 2 - ikC (k) [Ml] 


2b. 

i 


(A. 29) 


The [MO] and [Ml] matrices are defined in Eqs . (A, 11) and 
(A. 24). 

The distributed load resulting directly from the gust is given 
by Eq. (5o) : 


L 

_& - 




= 2 rTP a U 4 re£ a(y) |c(k) [J Q (k) - + iJ 1 (k)j 


= 27TP a Ub ref a(y)K(k) 
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The equivalent force is determined from 


£g}. 


27TP Ub l K 
a ref 


1 

(k) J (i - 


— j {q} ds 


= 2ttp Ub A K(k) a. 
a ref ' / i 


{GO} 


Ceil 


where : 


(A.?0) 



Assembling the matrices is straightf orward and is not performed 
here. The one boundary condition that enters in the assemblage pro- 
cedure comes from the assumption that the deflections are symmetric 
about the fuselage. This requires that the slope of the displacement 
be zero at the root, which in turn requires that the degree of freedom 
associated with this boundary condition be eliminated. 

A remaining task is the inclusion of the nons tructural masses of 
the fuselage and nacelle into the mass matrix. The fuselage is handled 
readily by approximating it as a point mass stationed at the root. This 


I 

l 


mass is simply added to the single element in the mass matrix that deals 
with the root displacement. 

The nacelle is slightly more complicated in that it is not posi- 
tioned at a node point of the finite element structure. The procedure 
used is to determine which element contains the nacelle and to then add 
to rhe mass matrix of that element terms corresponding to the equivalent 
mass of the nacelle obtained using a formula similar to that of Eq. 

1 

[m, nac " Vc f - W {nH ’' ,T ds 
0 

\ac ^( S NAC^ ^ s nac^ * (A. 31) 

where S( ) is the dirac delta. s„._ is the nondimertsionalized co- 
' ' NAC 

ordinate of the nacelle location. 
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